Single cell RNA sequencing after moderate traumatic brain injury: effects of therapeutic hypothermia

Abstract

Traumatic brain injury (TBI) initiates a cascade of cellular and molecular events that promote acute and long-term patterns of neuronal, glial, vascular, and synaptic vulnerability leading to lasting neurological deficits. These complex responses lead to patterns of programmed cell death, diffuse axonal injury, increased blood-brain barrier disruption, neuroinflammation, and reactive gliosis, each a potential target for therapeutic interventions. Posttraumatic therapeutic hypothermia (TH) has been reported to be highly protective after brain and spinal cord injury and studies have investigated molecular mechanisms underlying mild hypothermic protection while commonly assessing heterogenous cell populations. In this study we conducted single-cell RNA sequencing (scRNA-seq) on cerebral cortical tissues after experimental TBI followed by a period of normothermia or hypothermia to comprehensively assess multiple cell type-specific transcriptional responses. C57BL/6 mice underwent moderate controlled cortical impact (CCI) injury or sham surgery and then placed under sustained normothermia (37⁰C) or hypothermia (33⁰C) for 2 h. After 24 h, cortical tissues including peri-contused regions were processed for scRNA-seq. Unbiased clustering revealed cellular heterogeneity among glial and immune cells at this subacute posttraumatic time point. The analysis also revealed vascular and immune subtypes associated with neovascularization and debris clearance, respectively. Compared to normothermic conditions, TH treatment altered the abundance of specific cell subtypes and induced reactive astrocyte-specific modulation of neurotropic factor gene expression. In addition, an increase in the proportion of endothelial tip cells in the hypothermic TBI group was documented compared to normothermia. These data emphasize the importance of early temperature-sensitive glial and vascular cell processes in producing potentially neuroprotective downstream signaling cascades in a cell-type-dependent manner. The use of scRNA-seq to address cell type-specific mechanisms underlying therapeutic treatments provides a valuable resource for identifying targetable biological pathways for the development of neuroprotective and reparative interventions.

Introduction

Traumatic Brain Injury (TBI) is a major health concern and is one of the leading causes of death and disability worldwide [1, 2]. Every year an estimated 64–74 million people sustain a TBI [3]. Clinically, TBI can be diagnosed as mild, moderate, or severe based on the Glasgow Coma Scale (GCS) and considered a chronic health condition [4]. The pathogenesis of TBI is comprised of primary and secondary injury processes that have also been implicated in progressive neurodegenerative disorders such as chronic traumatic encephalopathy, Parkinson disease and Alzheimer’s disease [5,6,7,8,9]. While the primary brain injury directly leads to acute cell membrane disruption, vascular rupture, and blood-brain barrier (BBB) damage, secondary injury mechanisms include complex programmed cell death (PCD) and neuroinflammatory processes leading to delayed and progressive damage [10, 11]. Additionally, excitotoxicity, hypoxia, reactive molecular species, inflammation, and edema formation also contributes to progressive cellular damage [12]. In response to these multiple cell-specific injury process, the release of damage associated molecular patterns (DAMPs), cytokines and chemokines from damaged cells lead to the activation of microglia and astrocytes, and recruitment of circulating immune cells [13]. During this posttraumatic phase, the release of inflammatory mediators occurs from multiple cell types including neurons, microglia, endothelial cells (ECs), astrocytes, and peripheral immune cells [14]. Mechanical and other secondary injury mechanisms also alter endothelial structure and function leading to breakdown of the BBB and associated alterations in perivascular astrocytes and pericytes that participate in patterns of progressive injury [15]. In general, glial cell responses to neurotrauma include phenotypic changes indicating altered activation states assessed by morphological and genetic approaches as evaluated at both the regional and cell specific levels [16,17,18].

There have been advances in treatment of clinical TBI including acute surgical approaches to normalize intracranial pressure, reduce edema formation and secondary injuries. Also, pharmacological interventions directed at specific injury processes have been reported to promote various degrees of neuroprotection and improved behavioral outcome in experimental TBI models [19]. Unfortunately, no therapeutic strategies have been reported to improve neurological function in several randomized multicenter moderate/severe TBI trials [20, 21]. One strategy that has been translated to the clinic and evaluated in several neurological conditions such as TBI, SCI, stroke, cardiac arrest, and neonatal hypoxic encephalopathy is therapeutic hypothermia (TH) and targeted temperature management (TTM) [22,23,24,25,26,27]. While TTM is now commonly used to inhibit periods of hyperthermia in neuro intensive care units [28], there remains a lack of strong evidence to support the broad use of prophylactic hypothermia for improving neurological outcomes in moderate to severe TBI patients [11, 29,30,31,32].

Previous neurotrauma and cerebral ischemia studies have reported temperature-sensitive cellular and molecular changes associated with the benefits of TH in clinically relevant models [23, 33, 34]. These include the effects of postinjury brain cooling on necrotic and programmed cell death, BBB perturbations, inflammation, glial activation, cerebral blood flow and glucose metabolism [35,36,37,38,39,40]. Posttraumatic TH has also been reported to reduce molecular markers of cell death, glial activation and inflammatory process including the activation of the innate immune response to injury [35, 41,42,43]. These studies highlight the benefits of posttraumatic hypothermia on multiple cellular and molecular processes that may underlie the neuroprotective effects of TH after brain and spinal cord.

Recent advancements in single cell analysis including the spatial characterizations are clarifying cell-specific downstream signaling processes important for cell death, survival, and repair after CNS injury [44,45,46,47,48,49]. For example, several TBI studies have recently used genome-wide single cell RNA sequencing (scRNA-seq) analysis to describe cell type-specific changes and have characterized the phenotypic characteristics of subpopulations of specific cell groups [44, 46, 49,50,51]. In one study, Arneson and colleagues reported hippocampal cell-specific genes and signaling pathways 24 h after mild TBI and the effects of the T4 transporter transthyretin on these molecular changes [44]. In a follow-up study that conducted an in-depth analysis of the scRNA-seq data, Xing and colleagues identified several cellular and transcriptional responses to injury while emphasizing the importance of astrocytes and microglia in the acute inflammatory response to TBI [50]. In another study by Zheng and colleagues,13 different cell populations were identified with the upregulation of microglia and endothelial markers along with the downregulation of excitatory neurons reported [48]. Also, Witcher and colleagues focused on TBI-induced neuronal and inflammatory processes reporting that the majority of gene expression changes associated with inflammation and pathology were microglia independent at the subacute phase of TBI [52]. A recent study examined gene expression changes using 3 different models of TBI including a repetitive closed head model, controlled cortical impact (CCI), and a polytrauma CCI model [53]. This study found novel targets particularly for specific microglia and ependymal subtypes that vary based on TBI heterogeneity [53]. Taken together, these studies emphasize the power of scRNA-seq in characterizing cell specific responses and the opportunity to conduct hypothesis driven research to better understand which cell types and associated gene markers can serve as signatures of TBI pathology and therapeutic targets for therapeutic interventions. In this study we therefore compared the effects of early posttraumatic normothermia and hypothermia on cell-specific regional responses in several cell types. Specifically, we hypothesized that specific subpopulations of cells that participate in the reported benefits of cooling would be most sensitive to early posttraumatic hypothermia. To this end, we describe gene expression changes in neural, glial, and vascular cell clusters emphasizing the effects of temperature variations on microglial, astrocytic, and vascular and immune subtypes associated with neurotropic factor gene expression important in protective and reparative mechanisms.

Methods

Animals

All surgical procedures were conducted in accordance with the NIH guide for the care and use of laboratory animals and using animal protocols approved by the University of Miami Animal Care and Use Committee. Animals were housed in a temperature-controlled room (22⁰C) with a 12-h light/dark cycle. Animals were acclimated for at least 7 days prior to any experimentation. All animals had access to food and water ad libitum. We used 2- to 3-month-old male C57BL/6 mice, except one female Sham TBI + Hypothermia mouse.

Controlled cortical impact (CCI)

Animals were first weighed and then administered anesthesia via intraperitoneal injections (i.p) of ketamine (100 mg/kg) and xylazine (10 mg/kg). Animals then had their eyes covered with petroleum jelly to protect their eyes during the procedure. Hair was then removed from the animal’s skull utilizing Nair. Animals were next positioned into a stereotaxic frame and secured by ear bars and a bite plate. Once secured, a small incision was made using surgical scissors to expose the skull. A craniotomy was performed utilizing a surgical drill to drill a bore hole no deeper than the dura matter, bone fragments were cleared away using cotton tipped applicators. The craniotomy was located between lambda and bregma with a 4 mm diameter, 0.5 mm away from midline. Once the craniotomy was complete, the animal was placed under the piston of the CCI device. The piston was calibrated, and test fired prior to animal placement. CCI was performed between 9:00–10:00 am and hypothermia was given 30 min after injury. Once placed under the injury device, the animal was administered a moderate CCI (Custom Design and Fabrication, Richmond VA) utilizing a 4 m/sec velocity and 0.8 mm depth for 150 msec duration [54, 55]. At the conclusion of the injury, animals were removed from the stereotaxic device and excess blood was wiped away using cotton tipped applicators. Animals then had their incision closed with surgical staples and placed on a temperature-controlled pad for recovery. Animals were administered buprenorphine 0.1 mg/kg (s.c) post-operatively for pain management. Animals were observed during recovery and were returned to animal quarters upon awakening. Sham animals received all surgical interventions, respective drugs and temperature treatment but did not receive the craniotomy or TBI insult.

Temperature management

The following 4 groups were used for this study: TBI-Normothermia, Sham-Normothermia, TBI-Hypothermia, Sham-Hypothermia (n = 3 biological replicates per group; 12 replicates total; 1 animal/replicate). Therapeutic hypothermia was performed as previously described with some modifications (35). Briefly, systemic whole-body hypothermia was initiated 30 min after TBI or sham surgery to achieve target temperatures of 33^oC for hypothermia mice and 37^oC for normothermia mice. Cooling was initiated using body cooling packs with ice an once mice reached 33 ^oC rectal probes were used to monitor and sustain these temperatures for 2 h. Following temperature management, probes were removed, and mice were placed on heating pads with food and water for recovery after surgery.

Tissue dissociation for single-cell RNA-sequencing

At 24 hours post-injury, CCI or Sham animals were anesthetized with ketamine (100 mg/kg) and xylazine (10 mg/kg) and then transcardially perfused with oxygenated artificial cerebrospinal fluid (aCSF; 7 mM NaCl, 2.5 mM KCl, 1.25 mM NaH₂PO₄, 26 mM NaHCO₃, 75 mM sucrose, 20 mM glucose, 1.0 mM CaCl₂, and 2 mM MgSO₄). After extraction of the brain, the cortical injury site and surrounding penumbra was carefully dissected while avoiding underlying white matter tracts. Cortical tissue was chopped using a razor blade, washed with 5mL aCSF, and centrifuged at 300 g for 5min. After discarding supernatants, tissue pellets were enzymatically dissociated using the Miltenyi Neural Tissue Dissociation Kit-P (cat #130-092-628) according to manufacturer’s instruction with the following modifications. First, pellets were incubated in 2mL of enzyme mixture 1 for 30min at 37^oC (50µL Enzyme P in 1.9mL Buffer X) with gentle shaking by hand every 5min. 30µL of enzyme mixture 2 (10µL Enzyme A in 20µL Buffer Y) was added to the cell suspensions and manually triturated (slowly, 10x per sample) using a large-opening (1000µm diameter) fire-polished pipette. After 10min of incubation at 37^oC, suspensions were next triturated with a medium- and then small-opening (750µm and 500µm, respectively) fire-polished pipettes to produce single-cell suspensions. Suspensions were strained through a 70µm cell strainer, washed with 10mL aCSF, and centrifuged at 300 g for 5min. After discarding supernatants, pellets were incubated in 10µL Miltenyi Myelin Removal Beads II (Cat #130-096-733) and 90µL MACS buffer (0.5% Bovine Serum Albumin (BSA) in Hank’s Balanced Salt Solution (HBSS) without Ca²⁺/Mg²⁺) for 15min at 4^oC. Pellets were resuspended in 1mL MACS buffer and run through equilibrated LS Miltenyi MACS Magnetic Bead Columns (cat #130-042-401). Columns were washed with an additional 2mL MACS buffer and combined flow-through was centrifuged at 300 g for 5min. Pellets were resuspended in MACS buffer to achieve approximately 700–1200 cells/µL. Prior to library preparation, cell viability was assessed with ViaStain AOPI Staining Solution and only samples with viability > 75% were submitted. Approximately 10,000 cells/sample were loaded onto 10x Genomics Chromium and sample libraries were prepared using the Next GEM Single Cell 3’ Reagent Kit v3.1 (cat# PN-1000121) according to manufacturer’s instructions.

scRNA-seq data pre-processing and quality control

Sample libraries were sequenced on NovaSeq 6000 using paired-end sequencing with recommended cycle settings (Illumina, Inc.). The resulting raw base call files were demultiplexed into sample-specific FASTQ files and aligned to the Mouse mm10 (GENCODE vM23/Ensembl 98) reference genome using CellRanger v4.0.0 (10X Genomics). Further, we used the `cellranger count` command with default parameters to generate sample-specific feature-barcode matrices. Empty droplets were detected using a cell-calling algorithm with the `emptyDrops` function from the DropletUtils R package [56] to filter cell-containing droplets from the unfiltered `raw_feature_bc_matrix` outputs from CellRanger, which contain UMI counts for all 10x barcoded droplets. In brief, cells were ranked by total UMI count and visualized in a log-total UMI vs. log-rank plot. A spline curve was fit to the data to identify “knee” and inflection points. Cells with total UMI above the “knee” were considered cell-containing droplets while below the “knee” and above the curve inflection point were tested for statistically significant deviation in transcript profile compared to the ambient pool (total UMI < 200). To assess cell quality, we computed per-cell metrics such as total UMI, number of unique genes (defined as the total number of genes with UMI count greater than 0), and mitochondrial percentage (% of total UMI that maps to mitochondrial gene). A cell was defined to be of low-quality and removed from further analysis if it met any of the following criteria: (1) Its total UMI was less than 3 median-absolute-deviations below the per-sample median total UMI. (2) Its number of unique genes was less than an empirically determined threshold. For each sample, we estimated the probability density of the log10(# unique features) with the `density` base R function and calculated the slope of this probability density curve. We defined the minimum # unique gene threshold as the log10(# unique gene) value corresponding to the local minimum preceding the point of greatest slope of the log10(# unique gene) vs. slope curve. (3) Its mitochondrial percentage was greater than 25% i.e. more than 25% of the cell’s unique molecular identifiers (UMIs) mapped to a mitochondrial gene. After removal of low-quality cells, we proceeded to remove potential doublet cells (i.e. two or more cells captured in the same droplet). We applied the Scrublet Python package to each individual sample using default parameters. We also employed an empirical approach to remove doublets based on prior observations that algorithmic approaches could not reliably identify endothelial cell-pericyte doublets For each sample, we ran preliminary Seurat cluster analysis [57] using 20 principal components and FindClusters() `resolution = 1.3` and flagged putative doublet clusters based on well-known cell type-specific marker genes. FASTQ and associated files are available on Gene Expression Omnibus (GSE290150).

Normothermia scRNA-seq clustering analysis and annotation

For all clustering analyses, we used the Seurat Integration pipeline with modifications. After log normalization of gene counts, we used the `modelGeneVar` and `getTopHVGs` functions from the `scran` R package to identify genes varying within each sample the genes varying across all samples, respectively. In contrast to the Seurat default `FindVariableFeatures` which by default selects the top 2000 genes based on mean-variance modeling, the `scran` approach relies on statistical deviation of a gene’s variance from the mean-variance relationship to identify genes used for clustering analysis. After identification of variable genes, we performed sample-wise batch correction, data scaling, dimensional reduction, and clustering using Seurat functions. The number of principal components used as input for `FindNeighbors` and `RunUMAP` functions was determined using the elbow heuristic on a scree plot generated by `ElbowPlot`. Parameters for `FindClusters` were determined empirically by visualizing clustering stability with the clustree R package. Cluster data visualization was done with the UMAP (Uniform Manifold Approximation Projection) algorithm as implemented in Seurat. UMAPs and heatmaps were made with a color-blind friendly palette using scatterHatch [58].

For cell type/subtype annotations, we used a combined approach of differential expression (DE) testing, comparisons against reference data sets, and domain expertise of cell type-specific and subtype-specific marker genes. We ran DE testing with the `FindAllMarkers` Seurat function to identify marker genes and characterize clusters. We used default parameters, which implements a Wilcoxon rank-sum test comparing gene expression of cells within a given cluster versus all other cells. We repeated this test for each of the three cell compartments vascular, myeloid, neural) per treatment condition (TBI or hypothermia) to identify marker genes for each subtype. For comparisons against reference data sets, we used the SingleR algorithm to compare against mouse cortex scRNA-seq expression data by Zeisel et al. [59]. The `SingleR` method takes as input a reference dataset and a query dataset wherein the reference dataset has cells with annotated labels (e.g. by cell-type) and the query dataset has cells yet to be annotated. In brief, `SingleR` first performs pairwise DE tests to identify marker genes for each label in the reference dataset and then generates a prediction score for each query cell and each cell-type label based on the Spearman correlation of marker gene expression values between the reference and query data. To annotate major cell types, we used well-known markers to help annotate clusters corresponding to monocytes/macrophages (Ly6c2, Ms4a7, Ccr2), neutrophils (Mmp9, Ly6g), Microglia (P2ry12, Tmem119, Hexb), endothelial cells (Pecam1, Podxl, Cldn5), pericytes (Pdgfrb, Kcnj8, Cspg4), VSMCs (Acta2, Tagln), fibroblasts (Col1a1, Pdgfrb), choroid plexus cells (Ttr), ependymal cells (Foxj1, Dnah12), astrocytes (Gfap, Aqp4, Slc1a2), oligodendrocyte-lineage cells (Plp1, Olig2, Sox10), neural stem cells (Dcx, Ascl1), and neurons (Snap25, Thy1).

Prior to subtype cluster analysis, we ran Multi-Dimensional Scaling (MDS) analysis to quantify the contribution to gene expression heterogeneity due to categorical variables such as injury, temperature, cell-type, and source sample. This was done by computing the average expression profile per cell type by sample, computing a Euclidean distance matrix with the `dist` base R function, and generating an MDS plot with the `cmdscale` base R function.

For cell subtype cluster analysis, we performed a similar workflow as described above. For immune cell subtypes, we pulled and re-clustered monocytes/macrophages, neutrophils, and microglia. For vascular cell subtypes, we pulled and re-clustered endothelial cells, pericytes, VSMCs, and fibroblasts. For neural cell subtypes, we pulled and re-clustered choroid plexus cells, ependymal cells, astrocytes, OLCs, NSCs, and neurons. Prior to sub-clustering, we performed a first-pass clustering with Seurat’s standard workflow to identify any low-quality cells, doublets, or “contaminating” cell-types (e.g. Pecam1-expressing putative endothelial cells detected in the neural cell subtype analysis). Contaminant cell-types were identified using a list of well-described marker genes similar to the list presented above. We used the `FindMarkers` Seurat function to identify cell subtype marker genes as described above. UMAPs were created using `RunUMAP` as mentioned above. Neurons undergo cell death during the tissue dissociation process, so the remaining cells are highly skewed toward those that have survived tissue dissociation, and because of this bias, we excluded neurons from downstream analysis. However, we have included neuron cell numbers from each experimental group (Supp Fig. 1).

Fig. 1

Cellular composition of the cortical injury site after moderate TBI. (A) UMAP of cells from normothermic sham and TBI cortex at 24 h post-injury shows all major cell types (n = 3 biological replicates per group). Cells are colored by cell-type. (B) Same UMAP arranged by injury condition, sham (red) and TBI (yellow), shows TBI-specific cell types. (C) Proportion of each cell type before and after TBI shows cell type-specific effects of TBI, and that microglia, endothelial cells, and astrocytes comprise majority of the cell types in the dataset. Each data point is a biological replicate. Error bar = SEM (D) Dot plot of unique marker genes used to identify each major cell type

Analysis of effects of therapeutic hypothermia

To assess effects of hypothermia, we first used `SingleR` to classify cells from hypothermia-treated samples. Using annotated cells from the normothermia samples as reference dataset and unlabeled cells from hypothermia-treated samples as the query dataset, we applied `SingleR` to generate cell-type and cell subtype predictions. As a complementary approach, we also classified hypothermia-treated cells using Seurat’s `TransferData` function, which uses the “anchors” identified from the normothermia integrated analysis to generate predictions of cell-type labels in the query dataset. To assess the level of agreement between the two approaches for hypothermia-treated cells, we generated a multiple class confusion matrix with the Seurat-based prediction labels arbitrarily set as the “true” labels. To characterize the effects of hypothermia on cell gene expression, we identified differentially expressed genes (DEGs) using the `edgeR` R package according to recommended workflows [60]. For each cell-type, we generated pseudobulk expression data by summing counts grouped by sample and designed a separate `edgeR` model. We set the model matrix for the `glmQLfit` function to `~ (Injury + Temperature)^2` to capture the effects of “Injury” (i.e. TBI vs. Sham), “Temperature” (i.e. Normothermia vs. Hypothermia), and their interaction term (i.e. TBI: Hypothermia). DE test p-values were corrected using the Benjamin-Hochberg procedure by default. DEGs were visualized with volcano plots. For pathway analysis using DEGs, we used Qiagen’s Ingenuity Pathway Analysis (IPA) software. DEGs with a p-value ≤ 0.05 and an absolute logFC ≥ 1.0 were used as input. For Gene Ontology (GO) analysis, we used Fisher’s exact test as implemented in the topGO package. For comparison between two groups, we took all genes with an average log-fold change > 0.5 and adjusted p-value < 0.001.

In situ hybridization

Mice were anesthetized with avertin and transcardially perfused with cold 4% paraformaldehyde (PFA) in PBS. Brains were dissected, post-fixed in 4% PFA for two hours, and transferred to 30% sucrose in PBS solution for overnight incubation on a shaker. Injury sites (or corresponding regions in sham injured mice) were embedded in OCT compound (Tissue-Tek), and sectioned coronally at 16 μm thickness on a cryostat. Serial tissue sections were mounted onto Superfrost Plus slides and were stored at -20 °C. In situ hybridization was performed using the RNAscope Multiplex Fluorescent Kit v2 (Advanced Cell Diagnostics) according to the manufacturer’s protocol. Probes for Bdnf and Osmr were purchased from Advanced Cell Diagnostics. Following probe hybridization, each amplifier was individually incubated with slides for 30 min at 40 °C. After each individual fluorophore incubation, slides were washed with RNAscope 1X Wash Buffer and incubated with RNAscope HRP Blocker solution for 15 min at 40 °C. After in situ hybridization, slides were washed with PBS and incubated in DAPI for 5 min. Tissue sections were imaged on Andor Dragonfly spinning disk confocal microscope. ImageJ was used to count the number of Bdnf+/Osmr + nuclei around the injury site.

Results

scRNA-seq analysis of TBI reveals cellular heterogeneity of neural and non-neural cells

The early stages of moderate to severe TBI involve complex transcriptional changes among heterogeneous cell types. To gain a more comprehensive understanding of these processes, we used the CCI mouse model of TBI and performed scRNA-seq. C57BL/6 mice underwent either CCI or sham surgeries followed by either hypothermic or normothermic temperature management for 2 h using the following experimental groups: Sham Normothermia, Sham Hypothermia, TBI Normothermia, and TBI Hypothermia. At 24 h post-TBI, the cortical injury site and surrounding penumbra was carefully dissected for enzymatic and mechanical dissociation to obtain single cell suspensions which were then processed using the 10X Genomics Chromium platform. After full preprocessing and quality control, 88,200 cells across all 12 samples were used for downstream analysis (average 7350 cells per sample) with an average of 10,396 UMIs and 3023 unique genes per cell. An average of 10,397 UMIs were captured and an average of 3,023 unique genes were detected per cell (Supp Fig. 2).

Fig. 2

Vascular cellular heterogeneity after moderate TBI. (A) UMAP of vascular cells combined from sham and TBI normothermic animals (n = 3 biological replicates per group). Cells are colored by vascular subtype. (B) Same UMAP arranged by injury condition, sham (red) and TBI (yellow), shows TBI-specific cell types. (C) Proportion of each vascular cell type before and after TBI shows cell type-specific effects of TBI on endothelial cell subtypes (left) and perivascular cell subtypes (right). Each data point is a biological replicate. Error bar = SEM (D) Dot plot of unique marker genes used to identify each major cell type. (E) GO biological processes based on the top DEG of each subpopulation compared to all others and plotted on a log scale of its p-value (p < 0.001). The size of the circle represents the number of genes that defines the GO term, and the color of the circle represents its odds ratio

Previous research by Somebang et al. 2021 explored the diversity of immune cells in moderate to severe TBI [61]. However, to date, an in-depth analysis encompassing all cell types at the site of cortical injury has not been done. To fill this gap in knowledge, we concentrated our initial analysis on cells from normothermic samples − 23,606 cells from three normothermic TBI samples and 25,944 from three normothermic sham samples. We used Seurat’s CCA-based integration method to mitigate batch effects and subsequently performed a joint, unsupervised cluster analysis to identify shared and unique cell-types across injury conditions. Cells and clusters were visualized in a 2-dimensional projection with the UMAP algorithm (Fig. 1a). The cell-type identities of these clusters were determined by computing the top DEGs per cluster and cross-referencing them to previous scRNA-seq studies. We further reinforced these cell-type labels by applying the `SingleR` classification algorithm, with the Zeisel et al. 2015 mouse CNS dataset as a verified reference (Supp Fig. 3) [57]. This classification algorithm computes correlations between single cell expression profiles across datasets and allows for a data-driven cell-type labeling strategy, which showed good overall agreement with DEG-based labels. Our analysis revealed a total of 13 distinct clusters, encompassing all major neural and non-neural cell types typically found at the cortical site affected by CCI (Fig. 1a). These included neutrophils, microglia, various leukocytes, endothelial cells, pericytes, vascular smooth muscle cells (VSMCs), fibroblasts, choroid plexus cells, ependymal cells, astrocytes, neurons, neural stem cells (NSCs), and oligodendrocyte-like cells (OLCs). Notably, microglia, endothelial cells, and astrocytes, collectively comprising 70–80% of all captured cells, formed the most substantial clusters as shown in (Fig. 1c). Our differential expression analysis revealed novel, cell type-specific marker genes, such as Retnlg and Ly6c1, as well as marker genes that have been previously annotated (Fig. 1d). We also found certain subclusters of cells that were specific to the TBI injured group. These include a cluster of microglial cells, leukocytes, choroid plexus cells, ependymal cells, and neutrophils (Fig. 1a, b).

Fig. 3

Heterogeneity of neural cells after moderate TBI. (A) UMAP of neural cells combined from sham and TBI normothermic animals (n = 3 biological replicates per group). Cells are colored by neural cell subtype. (B) Same UMAP arranged by injury condition, sham (red) and TBI (yellow), shows TBI-specific cell types. (C) Dot plot of unique marker genes used to identify each major cell type. (D) Proportion of astrocyte subtypes before and after TBI shows cell type-specific effects of TBI. Each data point is a biological replicate. Error bar = SEM. (E) Volcano plot showing DEGs between GFAP-hi astrocytes (positive fold change) and reactive astrocytes (negative fold change). (F-G) GO biological processes based on the top DEG comparing reactive astrocytes (F) and Gfap-hi astrocytes (G) and plotted on a log scale of its p-value (p < 0.001). The size of the circle represents the number of genes that defines the GO term, and the color of the circle represents its odds ratio. (H) Ingenuity Pathway Analysis on the top DEGs between reactive astrocytes and sham astrocytes (Gfap-low and Gfap-hi). Pathways with p < 0.05 (threshold line on graph) shown

To better quantify the immediate impact of TBI on cellular diversity, we performed a multi-dimensional scaling (MDS) analysis, which attempts to decompose the variance in the data by categorical source such as cell-type, injury status, and sample source. The resulting MDS plots demonstrated that, while cell-type identity was the principal source of heterogeneity, TBI induced noticeable shifts in gene expression across multiple cell-types (Supp Fig. 3). Interestingly, the MDS plot also showed compartmentalization of cell-types according to vascular, immune, and neural ontogenies. From this, we reasoned that a sub-clustering analysis within each compartment would reveal additional cellular states and hypothesized that TBI induces unique cell subtypes across these compartments at 24 h post-injury.

Vascular cells after TBI

The vascular compartment was made of endothelial cells, fibroblasts, vascular smooth muscle cells (VSMCs), and pericytes. Taking these cells, we performed a sub-clustering analysis using a similar approach as before and visualized across injury conditions with UMAP (Fig. 2a) [45]. TBI caused changes in both gene expression and cell proportion in vascular cells. Among endothelial cells, we observed the appearance of tip cells, which was virtually absent in Sham mice, and a modest decrease in capillary ECs post-TBI, but no changes in arterial or venous ECs (Fig. 2b, c). Among the mural cells, we saw a reduction in pericyte proportion, the appearance of Fibroblast-1 and Fibroblast-2, and no change in VSMCs at this time point (Fig. 2c). Marker genes were then identified using differential expression analysis (Fig. 2d, suppl Fig. 4). Subtype analysis revealed a larger, aggregate cluster composed of Slc38a5(+) venous endothelial cells (ECs), capillary ECs with mixed expression, Sema3g (+) arterial ECs, and Apln (+) endothelial tip cells. Many of these markers are concordant with those identified in previous scRNA-seq studies analyzing vascular zonation in the brain [59] and the cellular response after spinal cord injury [45]. Our results also showed smaller clusters of pericytes, VSMCs, and two fibroblast subtypes. Pericytes were distinct from VSMCs, and fibroblasts based on expression of Kcnj8, a type of inwardly rectifying K⁺ channel. VSMCs were unique in their expression of Acta2, which encodes alpha smooth muscle actin. We also uncovered two fibroblast subtypes with unique expressions of Lum (encoding lumican, a sulfate proteoglycan) and Prg4 (proteoglycan 4).

Gene Ontology (GO) analysis of the top DEGs showed that arterial, capillary, and venous endothelial cells responded to TBI with an increase in genes associated with leukocyte migration, cell adhesion, and fibrinolysis, (Fig. 2e). Gene Ontology (GO) analysis of Tip Cell-enriched DEGs were associated with terms such as angiogenesis, cell migration, and actin filament polymerization. GO analysis of pericyte response to TBI revealed upregulation of genes related to hemostasis and regulation of angiogenesis and a decrease in aerobic respiration. Fibroblast-1 and Fibroblast-2 were both associated with collagen and extracellular matrix organization. Interestingly, direct comparison of the two subtypes showed that Fibroblast-1 was enriched for blood vessel development and adenylate cyclase signaling while Fibroblast-2 was enriched for metabolism and cellular respiration, suggesting heterogeneity in function among injury associated fibroblasts. Overall, our vascular cell analysis identified potential mediators of vascular remodeling as early as 24 h post-TBI.

Neural cells after TBI

Extraction and sub-clustering analysis of the neural compartment (choroid plexus cells, ependymal cells, astrocytes, OLCs, NSCs, and neurons, and excluding microglia discussed below) revealed 15 subtypes including various astrocyte types (GFAP-low, GFAP-high, reactive), radial glial-like cells, neural stem cells, neuroblast (NB), NB-derived neurons, oligodendrocyte progenitor cells (OPCs), OLCs, and specific neuron types each with their own distinct marker genes (Fig. 3a-c). This analysis revealed significant TBI-related changes in astrocyte subclusters (Fig. 3b, d). Heatmap for the top 20 DEGs for each neural cell subcluster is shown in supplemental Fig. 5. Prior research using mild TBI models showed acute phase astrocyte sensitivity [62]. Our study extends this to demonstrate single-cell level astrocyte changes in a moderate TBI model at 24 h post-injury. Notably, both GFAP-low and GFAP-high astrocytes, found in sham tissue, likely represent resting-state astrocytes. Post-injury, reactive astrocytes, expressing highest GFAP levels, emerged with a concomitant decrease in GFAP-low astrocytes compared to sham (Fig. 3d). Volcano plot analysis contrasting GFAP-high and reactive astrocytes identified several genes enriched in reactive astrocytes such as Hmox1, Serpina1a, and Osmr, with Osmr being highly specific to reactive astrocytes compared to all other neural cells (Fig. 3c, e). GO analysis of genes enriched in reactive astrocytes identified top biological processes related to cell proliferation and response to chemical stimulus, while the top biological processes in GFAP-high astrocytes were related to ATP synthesis and nucleoside triphosphate biosynthesis (Fig. 3f, g).

To further examine the mechanism of reactive astrocyte signaling after TBI, we performed Ingenuity Pathway Analysis (IPA) to identify pathways that were significantly upregulated in reactive astrocytes compared to Gfap-low and Gfap-hi astrocytes in the sham group (Fig. 3h). The most upregulated pathway was the Cytokine Storm Signaling Pathway, which indicates that astrocytes contribute to neuroinflammatory responses acutely after TBI. There was an upregulation of the Wound Healing and Acute Phase Response Signaling pathways, which is consistent with previously published studies showing the role of reactive astrocytes in wound contraction and limiting infiltration of immune cells after CNS injury [63]. We also observed in pathways such as Extracellular Matrix Organization and Integrin Cell Surface Interactions, which is consistent with the role of astrocytes in extracellular matrix deposition at the injury site [64, 65]. As expected, the pathways involving the neuroinflammatory response and cell death were upregulated, such as the interleukin-driven signaling circuits (IL-4, IL-13, and IL-6), macrophage-associated pathways, neuroinflammation signaling pathway, apoptosis signaling, and inflammasome signaling. On the other hand, upregulation of the Neurotrophin/TRK Signaling pathway within astrocytes after TBI shows their potential contribution to promoting neuronal survival and facilitating regeneration. These findings enhance our understanding of astrocyte-specific responses in the aftermath of TBI, shedding light on their complex roles and pinpointing potential targets for astrocyte-centered therapeutic intervention in TBI.

Myeloid cells after TBI

Circulating immune cell recruitment ensues almost immediately after TBI. As expected, TBI led to a significant increase in peripheral leukocytes such as neutrophils and monocyte/macrophages (Fig. 1c). MDS showed microglia were grouped closely with blood-derived leukocytes, speaking to their embryonic yolk sac ontogeny distinct from the neuroectodermal origin of other neural cells [66]. To further investigate this immune cell compartment, we extracted myeloid cells (monocytes/macrophages, neutrophils, and microglia), and performed sub-clustering, DE testing for marker genes, and UMAP visualization. Cluster analysis revealed subtypes of microglia as well as other myeloid cells including border associated macrophages, and dendritic cells (Fig. 4a, b). Additional heatmap analysis showing the expression levels of the top 20 immune cell markers in the Sham and TBI conditions is shown in supplemental Fig. 6. Most of the myeloid cells detected in our dataset were microglia, which separated into four subtypes (homeostatic, activated, disease-associated, and interferon microglia) each characterized by a distinct set of marker genes (Fig. 4a-c). Among microglia, the decrease in proportion of homeostatic microglia (higher expression of P2ry12, Gpr34, and Tmem119) after TBI was coupled with a concomitant increase in disease-associated microglia (lowest expression of the homeostatic genes and highest expression of Lpl, and Adam8), which were rarely detected in sham controls (Fig. 4d). However, both activated microglia and interferon microglia were detected in sham controls and proportions remained unchanged after TBI. The top 20 DEGs for each microglial subcluster is show in suppl Fig. 7. The presence of activated microglia in sham mice is most likely due to tissue dissociation-induced stress [67]. This is consistent with the likelihood that activated microglia represent an intermediate, transitional state between homeostatic and disease-associated microglia.

Fig. 4

Heterogeneity of myeloid cells after moderate TBI. (A) UMAP of neural cells combined from sham and TBI normothermic animals (n = 3 biological replicates per group). (B) Same UMAP arranged by injury condition, sham (red) and TBI (yellow), shows TBI-specific cell types. (C) Dot plot of unique marker genes used to identify each major cell type. (D) Proportion of myeloid subtypes before and after TBI shows cell type-specific effects of TBI. Each data point is a biological replicate. Error bar = SEM. (E) GO biological processes based on the top DEG comparing disease-associated microglia (left) and monocyte/macrophages (right) to all other myeloid cell subtypes and plotted on a log scale of its p-value (p < 0.001). The size of the circle represents the number of genes that defines the GO term, and the color of the circle represents its odds ratio

GO analysis of the top DEGs in disease-associated microglia revealed biological processes mostly associated with cell death and leukocyte migration (Fig. 4e). This could indicate that either disease-associate microglia are inducing cell death of nearby cells and/or are themselves undergoing cell death at this acute time point. GO analysis of monocytes/macrophages identified biological processes associated with leukocyte migration and glucose catabolism (Fig. 4e), which may indicate the migration of monocytes to the injury site and a metabolic shift to glycolysis during monocyte to macrophage differentiation.

Therapeutic hypothermia induces subtype-specific effects

Several studies have reported the effects of posttraumatic hypothermia on specific cell types including neurons, astrocytes, microglia, macrophages, oligodendrocytes, and cerebral endothelium [22, 23, 35, 36, 38,39,40,41, 43, 68,69,70]. However, no studies have reported the consequences of posttraumatic hypothermia using scRNA-seq to clarify the effects on specific cell subtypes and states at the single cell level. For this aim, we included six additional samples from TBI-hypothermia and Sham-hypothermia groups and performed cluster analysis. In our initial analyses, we found that the hypothermia-treated animals did not introduce any novel cell-types (Fig. 5a-c). Therefore, we implemented two orthogonal cell-type classification approaches– a transfer learning algorithm and a correlation-based algorithm. Using these methods allowed us to use the precise molecular markers we defined in our normothermia-only analysis to annotate hypothermia-treated cells. Overall, the two algorithms agreed strongly on cell subtype classifications. Given this strong concordance, we were next able to quantify precise changes in cell-type and subtype abundance after TBI and hypothermia and found that specific subtypes displayed altered proportions (Fig. 5d) and gene expression (Fig. 5e, suppl Fig. 9–11). TH led to reduced peripheral leukocytes including neutrophils and monocytes/macrophages, as well as neuroblasts, neural stem cells, ependymal cells, choroid plexus cells. Interestingly, TH resulted in an increased proportion of tip cells and reactive astrocytes. These different effects on subpopulation changes indicates that TH does not indiscriminately reduce cell proliferation or promote cell survival after TBI, but rather has cell type-specific effects.

Fig. 5

Effects of therapeutic hypothermia on cellular heterogeneity after TBI. A) UMAP of cells from normothermic and hypothermic sham and TBI cortex at 24 h post-injury shows all major cell types (n = 3 biological replicates per group (4 groups), total of 12 mice). B-C) Same UMAP of cells colored by major cell types (B) and treatment group (C) D) Proportion of cell subtypes in each of the four treatment groups. Each data point is a biological replicate. Error bar = SEM. E) Violin plot of the top five DEGs (by p value) between TBI + Normothermia and TBI + Hypothermia groups for each major cell type. Subtypes are in suppl Fig. 9–11

To assess the overall impact of TH at the transcriptional level, we compared the number of DEG between the major cell types and found that astrocytes were by far the most transcriptionally responsive to TH (Fig. 6b). Next, we compared the gene expression profiles of normothermic, and hypothermic reactive astrocyte subtypes visualized using a volcano plot (Fig. 6a). The comparison revealed many genes that were highly enriched in hypothermic reactive astrocytes, but Bdnf was highly specific to reactive astrocytes (and neurons) as shown in a feature plot of all cell types (Fig. 6c). Although previous study has shown that TBI itself can increase astrocyte Bdnf expression [71], our results suggest TH further increases Bdnf expression specifically in reactive astrocytes, and that this might be a key mechanism in TH’s neuroprotective role.

Fig. 6

Astrocytes are the most transcriptionally responsive to therapeutic hypothermia and upregulate BDNF. (A) Volcano plot of DEGs between hypothermia treated reactive astrocytes and normothermic reactive astrocytes. Top DEGs are labeled in red. (B) Quantification of the number of TBI-induced DEGs comparing normothermia and hypothermia for each major cell type. (C) UMAP feature plot of Bdnf expression across treatment groups. (D) Ingenuity Pathway Analysis of reactive astrocytes shows pathways that are affected by therapeutic hypothermia. Bars represent the -log(p-value) and color represents Z-scores; blue denotes negative Z-scores indicative of inhibition, and orange bars denote pathways with positive Z-scores indicative of activation. The threshold line indicates p < 0.05

IPA analysis of the top DEG in reactive astrocytes showed that the most highly affected pathways were downregulated after TH (Fig. 6d). The top pathway was glutathione-mediated detoxification followed by other well-known injury-related signaling pathways including HIF1-α, sphingosine-1-phosphate, eNOS, and Stat3. While many of these pathways have both beneficial and detrimental effects depending on the context and extent of activation, the overall downregulation of these pathways may reflect the astrocyte response to a less inflammatory (less infiltration immune cells) environment.

Validation of increased Bdnf in reactive astrocytes after therapeutic hypothermia

To further investigate the effect of TH on Bdnf expression in astrocytes, we used the single cell data to assess its expression across all major cell types and found that whereas Bdnf is only expressed in neurons before or after TBI, astrocytes were the only non-neuronal cell to express Bdnf after TH (Fig. 7a). To validate this in vivo, we used RNAscope to probe for Bdnf mRNA in tissue sections of the injury site at 1 day after TBI. To identify Bdnf expression in reactive astrocytes, we used Osmr as a second probe because this and prior studies have shown that it is virtually absent in homeostatic astrocytes but highly upregulated in reactive astrocytes [45], and because GFAP is expressed is non-reactive astrocytes. Although the single cell data shows that Osmr is also expressed in different vascular cell types (Fig. 7a), since Bdnf expression is specific to astrocytes, we compared Bdnf+/Osmr + nuclei between normothermic and TH groups. RNAscope showed visibly increased Bdnf signal around the injury site in TH brains. Although the cortical injury site was our focus, neuron-like cells in the hippocampus also showed increased Bdnf expression after TH (Fig. 7c). Quantification of Bdnf+/Osmr + nuclei in the perilesional areas showed a significant increase in TH mice compared to normothermic mice consistent with our scRNA-seq results (Fig. 7b).

Fig. 7

Validation of Bdnf expression using fluorescent in situ hybridization. (A) Dot plot of Bndf (left) and Osmr (right) expression in each major cell type in the four treatment groups. Osmr expression is used as a marker of reactive astrocytes. (B) Bar graph showing quantification of Bdnf + astrocytes (Bdnf+/Osmr+) in normothermic and hypothermic TBI cortex. n = 3 biological replicates per group. *p < 0.05, Student’s unpaired t-test. (C) RNAscope in situ hybridization to detect Osmr (red) and Bdnf (green) in hypothermic (left) and normothermic (right) TBI cortex. Cell nuclei labeled with DAPI (blue). Scale bar = 300 μm. (D) High magnification images of the injured cortex showing colabeling of Osmr and Bdnf. Scale bar = 15 μm

Discussion

TBI remains one of the largest causes of death and disability in the United States and results in over 300,000 hospital admissions and 2.8 million emergency department visits each year (CDC, 2014). Unfortunately, there are limited effective neuroprotective or restorative therapies available to mitigate severe post-TBI disability beyond surgical interventions or physical therapy. Over the past several years TH has been shown to reduce clinically relevant secondary injury mechanisms, thus providing neuroprotection and improved behavioral outcomes in models of TBI and other acute neurological conditions [19, 23, 72,73,74,75]. In terms of injury mechanisms underlying protection, studies using standard strategies for assessing morphological, biochemical, metabolic, and genomic changes have been utilized. For example, the use of TH has demonstrated reduction in indicators of posttraumatic inflammation using cDNA arrays and semiquantitative RT-PCR and documented families of genes that are influenced by this experimental therapy [33, 42]. Recent studies have used scRNA-seq as method for analyzing cell specific changes at the gene level in in various models of TBI and at both acute and chronic time points [52, 53]. Although there have been studies examining the cellular responses to posttraumatic TH, no temperature studies have yet evaluated patterns of single cell gene expression using scRNA-seq. We therefore examined the transcriptomic profile following either normothermic or hypothermic conditions to further understand cell type-specific changes that occur acutely after TBI and how treatment with TH might mitigate these changes.

In normothermic animals, TBI led to shifts in cell states that resulted in new subpopulations with distinct transcriptional identities at 24 h after injury. The most dramatic effects on cell state changes (i.e. virtually absent in shams, but present after TBI) were observed in astrocytes, microglia, and tip cells. The presence of activated microglia as well as Gfap-hi astrocytes in sham tissue could reflect the effects of dissociation-induced stress on these two glial populations. After TBI, there was a decrease in the Gfap-low subpopulation and an appearance of a reactive astrocyte subpopulation, which suggests that Gfap-low astrocytes gave rise to reactive astrocytes. However, graded changes in marker genes including Gfap, S1pr1, Btbd17, and Slo1c1 from Gfap-low to Gfap-hi to reactive astrocytes suggests that Gfap-hi astrocytes are a transition state between Gfap-low and reactive astrocytes. Similarly, TBI leads to decrease in homeostatic microglia and the appearance of a disease-associated microglia with activated microglia numbers unchanged after injury. Marker gene expression pattern showing a graded decrease in the canonical homeostatic microglia genes (P2ry12 and Tmem119) indicates that activated microglia are a transition state between homeostatic and disease-associated microglia, however, lineage tracing studies are necessary to properly test this hypothesis. Previous studies have used scRNA-seq to determine that microglia exhibit pro-inflammatory phenotypes at chronic time points [52]. Therefore, in future studies it will be important to examine chronic time points in TBI and TH model as well.

TH had different effects on different cell population after TBI. Consistent with previous studies, TH led to decreased presence of leukocytes including neutrophils and monocytes/macrophages [38, 41, 42, 76]. The most significant effect was on neutrophils, suggesting that limiting neutrophil infiltration may be a key inflammation target of TH. Our previous TBI studies demonstrated that TH for 4 h reduced microglial activation in a fluid percussion rat model at 24 h post injury [40]. However, our single cell data showed only a small decrease in the number of disease-associated microglia 24 h after TBI and TH suggesting potential species differences in the effect of TH on microglial activation states. An unexpected effect of TH was the increased percentage of tip cells and reactive astrocytes. The increase in these cell population versus a decrease in leukocytes indicates that the mechanism of action of TH is not via broad inhibition cell migration and/or proliferation but rather a more cell type-specific mechanism. There was a dramatic effect was on the increased percentage of tip cells. This could have been due to increased differentiation of endothelial cells to tip cells during angiogenesis, increased survival of tip cells, and/or reduced differentiation of tip cells to endothelial cells after TH. Similarly, the increased reactive astrocytes could have been due to increased transition from Gfap-low and/or Gfap-hi astrocytes, increased survival of reactive astrocytes and/or reduced reversion of reactive astrocytes to its resting state after TH. These possibilities remain to be tested in future experiments, and taken together provide additional insight into understanding the mechanisms of posttraumatic cooling and the complex cellular-dependent inflammatory responses after TBI.

It is interesting to note that some cell types, such as microglia and astrocytes, are more sensitive to tissue dissociation-induced stress than others, such as endothelial cells as indicated by the lack of activated endothelial subpopulations in the sham injured group. Although more time points are needed to better understand these difference effects, clues can be gleaned from the effects of injury on these cell types. Although TBI resulted in putative transition subpopulations in microglia and astrocytes, the only new endothelial subpopulation after TBI was tip cells without any evident pre-tip cell transition states. We made similar observations in our previous spinal cord injury scRNA-seq studies which included multiple acute time points [45]. Thus, we speculate that microglia and astrocytes are more sensitive to tissue dissociation-induced stress because they are more sensitive to CNS injury than compared to other cell types such as endothelial cells.

Since we loaded a target of 10,000 cells per sample, it does not accurately reflect the effect of TBI or TH on the cell number. In fact, although one would expect the cell numbers to decrease due to cell death after TBI, the cell numbers typically increase due to the presence of peripheral leukocytes that typically survive the dissociation protocol. This is why the percentage of each cell population is a better indicator than the number of cells. Nevertheless, our data shows that we had similar cell viabilities across treatment groups and samples, ranging between 76 and 87%, prior to loading (suppl Fig. 8) Furthermore, research suggests that males have a better neuroprotective response compared to females after brain injury and TH [77, 78]. Due to our low sample size, we had to limit our TBI groups to one sex. However further studies are needed at the single cell level to understand mechanisms for sex-specific TH management after TBI.

We found a significant increase in Bdnf in reactive astrocytes of TH treated mice after TBI. Although astrocytes have been reported as a major source of BDNF [79, 80], our data indicate very little BNDF expression by cells other than neurons both before and after TBI. This expression pattern is consistent with our previous scRNA-seq study on the normal and injured adult mouse spinal cord [45]. Interestingly, although TBI (or spinal cord injury) itself does not increase Bdnf expression, TH increases it significantly only in astrocytes. The mechanism underlying this astrocyte-specific effect of TH on Bdnf expression remains to be elucidated, but it could be a major mechanism of action of the benefits behind TH after CNS injury. BDNF delivered exogenously, or as a cell or gene therapy have shown beneficial effects after TBI in multiple animal models [81,82,83]. The astrocyte-specific increase in Bdnf after TH is likely to be more physiologically relevant and targeted than other therapeutic modalities. A recent study demonstrated that astrocytic BDNF plays a role in neurogenesis after CCI [84]. A future study examining astrocytic-specific manipulation of BDNF would enhance the understanding the mechanism of BDNF neuroprotection after TBI and TH. In our initial analysis, we did check for BDNF expression differences in all major cell types after TH and found a significant difference only in astrocytes. However, there does seem to be a trend in neurons in the sham groups, so it could be that our study is underpowered to detect this interesting effect.

Interestingly, we did notice an increase of BDNF expression in the hippocampal neurons after TBI and hypothermia, however very minimally in astrocytes. This supports previous studies demonstrating that TH increases BDNF expression in the ipsilateral cortex after TBI [85]. In the current study we examine cortical tissue at the single cell level, which is translationally relevant with our CCI model at the acute time point. In this study we only assessed the scRNA-seq changes in the injured cerebral cortex. Although we did not perform any quantifications, there were barely any BDNF + OSMR + astrocytes in the hippocampus in either the normothermic or the hypothermic animals. Future investigations assessing additional regional changes at various periods after TBI, both subacute and chronic, and TH treatments would be of interest to the field due to the reported effects of TH on brain regions remote from the primary contusive site. In addition, performing in situ hybridization or immunostaining for other cell types with BDNF and OSMR, such as neurons, would provide valuable information regarding cell-specific BDNF expression in various brain regions after TBI and TH.

Conclusion

In summary, we used single-cell RNA-seq analysis to compare the effects of posttraumatic normothermia and hypothermia on various cellular subpopulations that may be selectively sensitive to posttraumatic temperature modifications. This is the first study that tested the effects of posttraumatic temperature reduction on cell specific genetic responses that could participate in the reported neuroprotective effects reported in several injury models. We report that leukocytes, astrocytes and tip cells appear to be most responsive to TH in the subacute phase of injury. We also demonstrate that reactive astrocytes are significantly increased at this early posttraumatic time point, thus suggesting that TH may contribute to an increase in post-trauma tissue repair mechanisms after TBI. We have also identified novel pathways responsible for changes in astrocyte activity after TH treatment acutely after TBI. In addition, our data demonstrate that TH significantly increases the number of endothelial tip cells that also suggest that TH promotes angiogenesis acutely after TBI. Taken together, our transcriptomic results provide important information for the importance of understanding the temporal profile of the cell specific responses to TH treatment after TBI. This strategy may also show promise in elucidating individual cell groups and subpopulations that require attention when evaluating the potential success of experimental therapies that could show promise in the clinical arena. This finding is consistent with the differential effects of hypothermia on specific cell types and subpopulations in decreasing metabolic rates and conferring neuroprotection against secondary damage [86]. Collectively, these insights suggest that hypothermia may play a beneficial role in the TBI context by modulating astrocyte pathways that, while inherently part of the body’s response to injury, can lead to exacerbated damage or contribute to chronic neurological impairment if left unchecked.