Single cell molecular alterations reveal target cells and pathways of concussive brain injury
NATURE COMMUNICATIONS |
Single cell molecular alterations reveal target cells and pathways of concussive brain injury
Douglas Arneson 0 1
Guanglin Zhang 0
Zhe Ying 0
Yumei Zhuang 0
Hyae Ran Byun 0
In Sook Ahn 0
Fernando Gomez-Pinilla 0 2 3
0 Department of Integrative Biology and Physiology, University of California, Los Angeles , Los Angeles, CA 90095 , USA
1 Bioinformatics Interdepartmental Program, University of California, Los Angeles , Los Angeles, CA 90095 , USA
2 Department of Neurosurgery, University of California, Los Angeles , Los Angeles, CA 90095 , USA
3 Brain Injury Research Center, University of California, Los Angeles , Los Angeles, CA 90095 , USA
4 Institute for Quantitative and Computational Biosciences, University of California, Los Angeles , Los Angeles, CA 90095 , USA
5 Molecular Biology Institute, University of California , Los , USA
The complex neuropathology of traumatic brain injury (TBI) is difficult to dissect, given the convoluted cytoarchitecture of affected brain regions such as the hippocampus. Hippocampal dysfunction during TBI results in cognitive decline that may escalate to other neurological disorders, the molecular basis of which is hidden in the genomic programs of individual cells. Using the unbiased single cell sequencing method Drop-seq, we report that concussive TBI affects previously undefined cell populations, in addition to classical hippocampal cell types. TBI also impacts cell type-specific genes and pathways and alters gene co-expression across cell types, suggesting hidden pathogenic mechanisms and therapeutic target pathways. Modulating the thyroid hormone pathway as informed by the T4 transporter transthyretin Ttr mitigates TBI-associated genomic and behavioral abnormalities. Thus, single cell genomics provides unique information about how TBI impacts diverse hippocampal cell types, adding new insights into the pathogenic pathways amenable to therapeutics in TBI and related disorders.
T sports, and military environments and often leads to
longraumatic brain injury (TBI) is common in domestic,
term neurological and psychiatric disorders1. The
hippocampus is a member of the limbic system and plays a major role
in learning and memory storage. As a major aspect of the TBI
pathology2, hippocampal dysfunction leads to memory loss and
cognitive impairment. The hippocampal formation encompasses
four Cornu Ammonis (CA) subfields largely composed of
pyramidal cells, and their connections with dentate gyrus (DG) cells.
The CA?DG circuitry has served as a model to study synaptic
plasticity underlying learning and memory. Glial cells are vital to
the hippocampal cytoarchitecture, however, their interactions
with neuronal cells are poorly defined. The heterogeneous
properties of the hippocampal cytoarchitecture have limited the
understanding of the mechanisms involved in the TBI pathology.
Mild TBI (mTBI) is particularly difficult to diagnose given its
broad pathology, such that there are no accepted biomarkers for
mTBI3. This limitation becomes an even more pressing issue
given the accumulating clinical evidence that mTBI poses a
significant risk for neurological and psychiatric disorders associated
with the hippocampus such as Alzheimer?s disease (AD), chronic
traumatic encephalopathy (CTE), post-traumatic stress disorder
(PTSD), epilepsy, and dementia4. Accordingly, there is an urgent
need to identify functional landmarks with predictive power
within the hippocampus to address current demands in clinical
Given that gene regulatory programs determine cellular
functions, scrutiny of large-scale genomic changes can reveal clues to
the molecular determinants of mTBI pathogenesis including
cellular dysfunction, injury recovery, treatment response, and
disease predisposition. However, existing genomic profiling
studies of mTBI are based on heterogeneous mixtures of cell
conglomerates5?9 which mask crucial signals from the most
vulnerable cell types. Here, we report the results of a high
throughput parallel single cell sequencing study, using Drop-seq,
to capture mTBI-induced alterations in gene regulation in
thousands of individual hippocampal cells in an unbiased manner. We
focus on concussive injury, the most common form of mTBI,
using a mild fluid percussion injury (FPI) mouse model which
induces identifiable hippocampal-dependent behavioral deficits
despite minimal cell death10. We examine the hippocampus at 24
h post-mTBI, as this is a pivotal timeframe for pathogenesis and
is generally used for diagnostic and prognostic biomarker
To our knowledge, this is the first single cell sequencing study
to investigate the mTBI pathogenesis in thousands of individual
brain cells in parallel, offering a cell atlas of the hippocampus
under both physiological and pathological conditions. In doing
so, we provide novel evidence about the cellular and molecular
remodeling in the hippocampus at the acute phase of TBI and
help answer critical longstanding questions. Which cell types are
vulnerable to mTBI at the acute phase? Within each cell type,
which genes have altered transcriptional activities that are
induced by mTBI? Which molecular pathways are perturbed by
mTBI in each cell type and how do they relate to mTBI pathology
and pathogenesis of secondary brain disorders such as AD and
PTSD? How do the coexpression patterns of genes across cells
and circuits vary in response to mTBI? Through answering these
questions, the identified sensitive cell types and associated gene
markers can serve as signatures of mTBI pathology that inform
on the stage, functional alterations, and potential clinical
outcomes. Since the cell is the elementary unit of biological structure
and function, we reveal fundamental information that can lead to
a better understanding of the mechanistic driving forces for mTBI
pathogenesis and identify potential pathways for therapeutic
targeting in an unbiased manner. As a proof of concept, we use
the data-driven single cell information to prioritize Ttr, encoding
transthyretin, and show for the first time that modulating Ttr
improves behavioral phenotypes and reverses the molecular
changes observed in mTBI.
Unbiased identification of cell identities in hippocampus.
Using Drop-seq12, we sequenced 2818 and 3414 hippocampal
cells from mTBI and Sham animals, respectively. A single-cell
digital gene expression matrix was generated using Drop-seq
Tools12 and subsequently projected onto two dimensions using
tdistributed stochastic neighbor embedding (t-SNE)13 to define
cell clusters (Methods). We detected 15 clusters each containing
cells sharing similar gene expression patterns (Fig. 1a). The cell
clusters were not due to technical or batch effects (Supplementary
To resolve the cell-type identities, we obtained cluster-specific
gene signatures (Supplementary Data 1) and compared them to
known signatures of hippocampal cell types derived from
Fluidigm-based single cell studies14,15 (Methods). We recovered
all known major cell types including neurons, oligodendrocytes,
oligodendrocyte progenitor cells (oligoPCs), microglia, mural
cells, endothelial cells, astrocytes, and ependymal cells (Fig. 1b).
Previously known cell markers, such as Aqp4 for astrocytes, Mog
for oligodendrocytes, and C1qc for microglia, all showed distinct
cluster-specific expression patterns, confirming the reliability of
our data-driven approach in distinguishing cell types (Fig. 1c?e;
additional examples in Supplementary Fig. 2). Beyond retrieving
known cell markers, we also identified novel marker genes for
each hippocampal cell type, such as Calml4 for ependymal, Vcan
for oligodendrocytes, and Ly6a for endothelial cells (Fig. 1f;
Supplementary Data 1). We confirmed the cellular identity and
hippocampal localization of many novel marker genes using the
in situ hybridization (ISH) data from the Allen Brain Atlas (Fig. 2;
Supplementary Figs. 3?9).
We identified two potential novel cell clusters whose gene
expression signatures did not significantly overlap with any of the
known hippocampal cell types, and were hence named
Unknown1 (113 cells; 1.8% of all cells tested) and Unknown2
(74 cells; 0.9%; Fig. 1a). Unknown1 has markers indicative of cell
growth and migration such as Ndnf, Nhlh2, Reln, and Igfbpl1, as
well as endothelial markers (Fig. 1b), suggesting the cells are
migrating endothelial cells. Endothelial cells as main components
of the blood brain barrier which is disrupted after mTBI16, with
proliferation and migration of endothelial cells being an intrinsic
aspect of new vessel formation17. Unknown2 expresses unique
markers indicative of cell differentiation such as Pcolce, Col1a2,
Asgr1, Serping1, and Igf2, along with markers of endothelial,
mural, and ependymal cells (Fig. 1b), suggesting that these may
be progenitor cells differentiating into multiple lineages. To
further resolve the identities of these unknown clusters, we
examined the expression patterns of their signature genes in the
Allen Brain Atlas ISH data. The marker genes of the Unknown2
cluster colocalized to a population of cells in the choroid plexus
distinct from the ependymal cells (Fig. 2; Supplementary Fig. 9).
The localization of these cells to an area implicated to house stem
cells and progenitor cells18 in conjunction with the functions of
the marker genes, supports that this cluster may represent
progenitor cells. These cell clusters may illustrate cell types that
have been previously left undetected by classical morphological
categorization, although detailed functional characterization is
required to confirm their identities and functions.
To further characterize neuronal diversity within the
hippocampus, we took the cell cluster that expresses canonical neuronal
markers (Fig. 1a) and further refined nine subclusters using the
CA Subtype1 neurons
CA Subtype2 neurons Ependymal Astrocytes CA1 neurons
BackSPIN biclustering method which offers better resolution14
(Fig. 3a). Annotation with known neuronal markers helped
resolve GABAergic interneurons, DG granule cells, and 4 subtypes
of CA pyramidal neurons (Fig. 3b). However, two clusters
(Neuronal Subtype1, Neuronal Subtype 2) express general
neuronal markers, but do not express markers of any currently
established neuronal subtype in the hippocampus. The functions
of their markers (Fig. 3b) suggest that these clusters may contain
neurons with the potential to differentiate or self-renew. The
unique ability of Drop-seq to catalog cells based on unbiased
genomic parameters was also instrumental in unveiling novel
neuronal markers14,15 including Ptn in CA1 neurons, Ly6e in
CA3 neurons, and Sema5a in DG granule cells (Fig. 3b; full list in
Supplementary Data 1). We confirmed the subtype-specific
expression patterns in hippocampal subregions and cell types
(Fig. 3c?e) and further verified the expression specificity of
several novel markers using the Allen Brain Atlas ISH data (Fig. 2,
Supplementary Figs. 3?9). The CA Subtype2 cells were found to
be located in the Subicular Complex (Fig. 2), which mediates the
main output of signals from the hippocampus. Like the CA
subregions, the Subicular Complex is also comprised of pyramidal
neurons, which explains the expression of CA neuronal signature
genes in addition to the genes specific to this cell cluster. The
subiculum is involved in pathology of CTE and the associated
dementia following TBI4.
These results indicate that our transcriptome-driven, unbiased
Drop-seq approach has the unique ability to uncover potential
new cell types, states, and markers based on genomic features that
may infer function. Based on the central dogma, gene regulation
is upstream of the production of proteins, which are fundamental
for cell function. In contrast, morphology-based approaches may
not offer the resolution to distinguish subtypes of cell populations
that share similar morphology but carry certain unique functions.
For instance, cells of the two CA subtype clusters uniquely
express specific marker genes which may infer unique functions.
However, detailed functional characterization of these potential
new cell subtypes is required in future studies to test functional
differences and their role in the hippocampal circuitry under
homeostatic and/or pathological conditions.
Identification of vulnerable hippocampal cell types to mTBI.
To retrieve hippocampal cell types vulnerable to mTBI, we first
compared the cell population proportions between mTBI and
Sham animals. Ependymal cells were found to be more abundant
in mTBI compared to Sham animals (89% from mTBI samples vs
11% from Sham samples, Fig. 1a). This large shift in relative
abundance of ependymal cells at 24 h post-surgery is consistent
with the reported role of ependymal cells in acute post-injury
processes such as circuit repair, scar formation, and potential
source of progenitor cells19,20. We also see significant increases in
the cell proportions of microglia, endothelial cells, and Unknown
2, as well as significant decreases in neurons, oligodendrocytes
and oligoPCs post-TBI (Supplementary Table 1). However, many
experimental factors in the Drop-seq procedure may influence the
capture rate of different cell types between samples, and changes
in the relative proportion of a cell type do not directly implicate
cell proliferation or death but could be the result of shifts in other
cell types. Therefore, caution is needed in the interpretation of
We also examined the shifts in global transcriptome patterns
within each cell type by testing if the Euclidean distance of gene
expression profiles of cells between the two groups is larger than
expected by chance (details in Methods). This analysis revealed
that DG, ependymal cells, astrocytes, microglia, oligodendrocytes,
oligoPCs, endothelial cells, Neuronal Subtype 1, CA1, CA
Subtype1, and GABAergic neurons demonstrate significant global
transcriptomic shifts between mTBI and Sham (Supplementary
Fig. 10). In particular, mTBI had such a profound effect on the
transcriptome of DG granule cells, that they became clearly
separated into two distinct clusters (Figs. 1a, 3a). Post-traumatic
epilepsy is a major concern in TBI and is attributable to DG
The above analyses implicated the majority of hippocampal cell
types to be vulnerable to mTBI and open the possibility for
investigating the specific roles of each cell type in mTBI pathology
and for designing targeted treatments. For example, the genomic
markers of the DG granule cells highly sensitive to mTBI can be
used to design treatments that target specific DG subpopulations
responsible for post-traumatic epilepsy to avoid the side effects of
classical epilepsy treatments targeting broad populations of cells.
mTBI alters cell-cell gene co-expression in the hippocampus.
Emerging evidence in the neuroimaging field suggests that
changes in the interaction patterns among cells in circuits can
contribute to reduced cognitive capacity in TBI22. We used the
co-regulation patterns between genes of different cell types to
infer interactions among cells, as gene co-expression can infer
functional connectivity23,24. Specifically, we focused on marker
genes encoding secreted peptides from each cell type (source
cells), which have the potential to interact with receptors and
regulate downstream effector genes in other cell types (target
cells) (Methods; Fig. 4a). This approach was found to recapitulate
known cell?cell communications within the trisynaptic circuit of
the hippocampus (Fig. 4b?d). The gene co-expression patterns
reflected the known mutual interactions between DG and CA3 as
well as the CA3 to CA1 interaction. We applied our method to
the mTBI data and found extensive reorganization in the pattern
of gene coexpression, potentially representing interactions,
among cells in response to mTBI. For example, the interaction
from astrocytes and ependymal cells to neurons and from
microglia to oligodendrocytes was enhanced in mTBI (Fig. 4e, f).
We also found decreased interaction between microglia and
neurons, and decreased interaction from oligodendrocytes to
neurons in mTBI. These shifts in the pattern of gene
coregulation among hippocampal cell types may implicate
reorganization of the working flow in response to mTBI challenge. The
reorganization was also reflected in changes in the association
pattern among single genes. For instance, the correlation of Bdnf
from neurons with cell metabolism genes in microglia was lost in
mTBI. Notably, the known AD risk gene Apoe in astrocytes and
ependymal cells (source cells) showed strong correlations with
mitochondrial metabolism genes in neurons (target cells) after
mTBI. These results suggest the potential of mTBI to alter the
interactive patterns at the level of circuits and genes in the
hippocampus, with cell metabolism involved in the interactions
(detailed peptide-gene correlations in Supplementary Data 2). We
acknowledge that this in silico analysis points towards potential
circuit perturbations but would require downstream validation.
Stxbp6 Nts Odf3b Mgp
TBI DG granule cells
Neuronal subtype 2
Sham DG granule cells
CA1 neurons: Wfs1
DG granule cells: Dsp
CA subtype 2 neurons
GABAergic interneurons: Gad2
CA subtype 1 neurons
Neuronal subtype 1
CA1 neurons (139)
CA3 neurons (64)
CA subtype 1 neurons (67)
CA subtype 2 neurons (77)
Sham DG granule cells (199)
TBI DG granule cells (218)
GABAergic interneurons (223)
Neuronal subtype 1 (78)
Neuronal subtype 2 (
Corrected P value
CA subtype 1 neurons
CA subtype 2 neurons
DG granule cells
Neuronal subtype 1
Neuronal subtype 2
uniquely points to the specific cell types engaging these pathways
and offers novel insights into the functions of individual cell
types, including previously understudied cell populations, in
TBI is followed by a stage of metabolic dysfunction28 which
reduces the capacity of the brain for activity-dependent
plasticity29. Our results provide a better understanding of the
molecular and cellular underpinnings of the metabolic crisis
typical of mTBI by identifying down-regulation of mitochondrial
metabolic genes in astrocytes and CA1 pyramidal cells. Astrocytes
supply energy substrates to neurons and are essential for neuronal
function. Interestingly, our results indicate that CA1 pyramidal
cells also experience metabolic alterations at this early stage of
mTBI pathogenesis (24 h). In addition, DEGs in CA1 pyramidal
cells informed on increased expression of glutamate transporters,
which could explain the altered capacity to sustain long-term
The risks posed by mTBI on the development of other
neurological disorders are a pressing issue in clinical
neuroscience. The genes, pathways, and the associated cell types that
we identified provide insights into the molecular and cellular
bases for these disease connections. DG granule cells, which
function to interact with CA pyramidal cells (Fig. 4b), showed
alterations in genes involved in cell-cell signaling (Npy, Penk,
Ptprn, and Ihnba) as well as neuroplasticity genes such as Bdnf
and Ntf3. These pathways may underlie the sensitivity of DG
granule cells to contribute to post-traumatic epilepsy after mTBI.
Pathways informed by DEGs from the endothelial and ependymal
cells implicate the importance of these cell types in amyloid
buildup, as indicated by the upregulation of the known
proamyloid deposition gene B2m and down-regulation of inhibitors
of beta-amyloid aggregation (Apoe, Itm2a, Itm2b, and Itm2c) in
these cell types. As dysregulated metabolism and amyloid
deposition are key features in AD, CTE, and PD31, our study
provides detailed information on the specific cell types such as
astrocytes, CA1 pyramidal cells, endothelial and ependymal cells
that could be the starting loci for the wave of post-mTBI amyloid
buildup in chronic mTBI.
mTBI alters gene expression in a cell-type specific manner. We
found that many of the DEGs were significantly altered by mTBI
in a specific cell type (Fig. 5a, b) and showed clear cell
typespecific shifts in expression patterns (Fig. 5c?e). Importantly,
>50% of the DEGs identified at the single cell level were masked
in bulk tissue analysis (Fig. 5f; tissue-level DEGs in
Supplementary Data 3) and these unique cell-level DEGs were primarily
from cell types of low abundance such as neuronal
subpopulations. On the other hand, the common DEGs between single-cell
and tissue-level analyses were mainly from abundant cell types
such as astrocytes and oligodendrocytes. As less abundant cells
such as neurons carry essential functions in the brain, the cell
type-specific DEGs can be strong drivers of disease but will be
missed in bulk tissue analysis. Therefore, genomic information in
individual cell types has the advantage to extract hidden
mechanisms involved in TBI pathology that otherwise would be
masked in bulk tissue studies.
Cell type-specific DEGs may serve as selective biomarkers or
therapeutic targets that can trace or normalize specific
abnormalities of mTBI pathology (Fig. 6a). For instance, Id2, a gene
previously described to be upregulated by seizures in DG32, is
specifically upregulated by mTBI in DG granule cells and could
serve as a potential target to temper post-traumatic epilepsy. We
also found that P2ry12, a marker for brain resident microglia33, is
specifically downregulated in microglia post-mTBI, suggesting the
potential utility of P2ry12 as a marker for early inflammatory
response to TBI. Interestingly, Trf, encoding transferrin for iron
transport was upregulated in oligodendrocytes, suggesting a
possible involvement of Trf on the described association between
iron deposition and cognitive deficits in mTBI34,35. Several of the
cell-type specific DEGs are related to amyloid deposition and AD,
including Apoe?an ependymal-specific DEG and known for its
effects on AD and TBI36, and Itm2a?an endothelial specific DEG
and an inhibitor of amyloid-beta production and deposition37.
These results indicate that the putative action of TBI in amyloid
regulation involves various cell types.
We also found several cell-type specific mTBI responsive genes
that have not been implicated in mTBI previously. For instance, a
CA1 pyramidal cell-specific DEG Klhl2, which encodes an actin
binding protein, was recently implicated in human neuroticism38,39
and showed increased neuronal expression post-mTBI (Fig. 6a).
Klhl2 in CA1 pyramidal cells may serve as a novel link between
mTBI and the increased tendency for neuroticism observed in blast
TBI40. Arhgap32 is primarily upregulated in CA subtype 2 cells in
mTBI. It encodes a neuron-associated GTPase-activating protein
that may regulate dendritic spine morphology and strength41.
Arhgap32 may be critical for supporting transmission of
information across hippocampal cells. The endothelial-specific gene Fxyd5
encodes a glycoprotein that enhances chemokine production and
inhibits cell adhesion by downregulating E-cadherin42. It was
upregulated in mTBI in endothelial cells, suggesting that it may play
a role in neuroinflammation and blood-brain-barrier dysfunction
associated with TBI.
DEGs listed as up had increased expression in TBI and DEGs listed as down had decreased expression in TBI compared to sham. Cell types shown in parenthesis are rare cell types with fewer cells
analyzed, and hence did not have sufficient numbers of DEGs at FDR < 5% between Sham and TBI samples to reveal significant pathways. Instead, DEGs reaching a threshold of p < 0.01 (number of DEGs
shown in parenthesis) were used to derive suggestive pathways for these rare cell types. P values were determined using a bimodal likelihood ratio test
Validation of cell-type specific mTBI DEGs. We used the
RNAscope multiplex RNA ISH to validate select cell-type specific
mTBI DEGs (selected from Fig. 6a) by co-hybridizing
hippocampal tissue slices with a cell identity marker and a DEG for
each select cell type. Using n = 8 animals per group, over 900
zstack images with two ISH channels (one for cell markers and the
other for DEGs) and a DAPI (for nuclei) channel were generated
for multiple DEGs and cell types. RNAscope is a high-resolution
quantitative ISH which detects single mRNA molecules.
Colocalization between a cell marker gene and a DEG is determined
by the presence of fluorescent spots representing single mRNA
molecules of both genes within the same cell. Automated cell
segmentations were performed using the CellProfiler software43
and quantification of mRNA molecules per cell was carried out
using FISH-quant44. To quantify cell type-specific differential
expression of a DEG, we first determined the cell type of interest
from a z-stack image by thresholding the counts of the cell type
marker, followed by comparing the counts of the DEG mRNA
molecules only within that cell type between Sham and TBI
samples (Methods). As shown in Fig. 7 (high magnification
images of select regions) and Supplementary Fig. 11 (low
magnification images of larger areas), the select DEGs are localized in
the same cell segments expressing the corresponding cell identity
markers. Quantitative assessments validated significant
differential expression of cell type-specific DEGs in Microglia (P2ry12),
DG granule cells (Id2), and CA Subtype 2 neurons (Arhgap32) in
both CA1 and CA3 regions (p < 0.05 by bimodal likelihood ratio
test; n = 8/group). Klhl2 in CA1 pyramidal neurons did not reach
statistical significance (p = 0.07) but showed an increasing trend
which is consistent with the Drop-seq data.
Prioritization of mTBI targets based on single cell data. The
aforementioned cell-type specific genes perturbed by mTBI
provide unique information about the microcircuits underlying the
mTBI pathophysiology. These can be leveraged for the design of
therapeutic interventions to target specific cell types driving
pathological manifestations if their causal roles in pathogenesis
are confirmed. Conversely, identifying DEGs that are affected
across multiple cell types by mTBI has the potential to pinpoint
the most vulnerable genes that are responsible for the broad
symptomology of mTBI. Such genes cannot be retrieved without
examining multiple cell types individually to confirm the
widespread effect across cell types. Given the broad aspects of the TBI
pathogenesis, targeting such pan-hippocampal DEGs may offer
broader and stronger therapeutic effects by normalizing the
functions of multiple cell types.
Based on the consistency of differential expression across cell
types, we prioritized a number of multi-cell type DEGs affected
by mTBI (Fig. 6b). Several such DEGs are related to beta-amyloid
and AD, including Ttr?encoding transthyretin which is an
amyloid beta scavenger45, mt-Rnr2?encoding the mitochondria
factor humanin46 which is protective against AD, Itm2b?an
inhibitor of beta-amyloid deposition47, and Apba2 (Mint2)?a
stabilizer of amyloid precursor protein APP48. These
panhippocampal DEGs, along with some of the cell-type specific
DEGs previously discussed, strongly implicate pathways related
to post-mTBI AD pathogenesis and can be targeted for
postmTBI AD prevention. Slc17a7 is a pan-neuronal DEG. It encodes
a sodium-dependent phosphate transporter in neuron-rich
regions of the brain and functions in glutamate transport49.
Interestingly, genetic polymorphisms in this gene have been
associated with recovery time and severity of outcomes after
sport-related concussion in humans50.
Notably, the gene Ttr represents the most robust DEG across
hippocampal cell types in that it was a top DEG with increased
expression post-mTBI in 7 of the 10 major cell types and 6 of the
8 neuronal clusters, with the highest expression and strongest
induction seen in ependymal cells (Fig. 6b). Therefore, the full
information across hippocampal cell types established Ttr as the
most consistent DEG post-mTBI, providing a strong rationale to
82 74 71
prioritize Ttr for testing. Ttr encodes transthyretin, a transporter
that carries the thyroid hormone thyroxine, preferentially T4,
across the blood brain barrier51. We hypothesize that the
panhippocampal upregulation of Ttr might implicate an impaired
thyroid hormone pathway in TBI and a compensatory need for
thyroxine T4, the major brain-specific substrate of transthyretin.
Given the strong implication of altered cell metabolism in various
cell types discussed above, the critical role of thyroid hormone in
regulating metabolism could serve as a platform to mitigate the
metabolic crisis after mTBI across a broad range of hippocampal
As a proof of concept, we examined the possibility that
modulating the thyroid hormone pathway and Ttr may serve as
an effective strategy to counteract mTBI pathology. We first
quantitatively confirmed, using RNAscope multiplex ISH assays,
that Ttr was indeed significantly upregulated in multiple cell types
and brain regions as predicted by Drop-seq (Fig. 8 for high
magnification images, Supplementary Fig. 12 for low
magnification images). Acute intraperitoneal injection of T4 within the first
6 h post-mTBI protected learning and memory, as determined by
the Barnes Maze test one-week post-mTBI. For the learning
component, animals were trained with two trials per day for four
consecutive days, and memory retention was assessed two days
after the last learning trial. Differences in learning and memory
between groups were determined using a two-tailed Student?s
ttest. The learning effects were evidenced by a sustained reduction
in the latency to find the escape hole across all time points in the
T4 group compared to TBI mice without T4 (Fig. 9a, b). The
effects of T4 on memory were demonstrated by the recovery of
the latency time to a level comparable to the Sham group (Fig. 9c,
d). The improved Barnes Maze performance was not due to
changes in velocity (Supplementary Fig. 13). It is noteworthy that
T4 was injected immediately after the lesion, which could protect
a variety of cellular functions against the disruptive effects of TBI.
This is the first time that T4 treatment was found to mitigate
cognitive deficits in a mouse model of concussive injury.
A recent study using cortical impact injury suggested a
potentially protective action of T4 by testing candidate genes
related to hypoxia and neurogenesis52. Our analysis of known T4
transporters indicates that T4 treatment primarily downregulates
Ttr and has weaker effects on genes encoding other transporters
(Fig. 9e). Genes encoding thyroid hormone receptors that are
downstream of T4 also had less consistent changes across cell
types compared to Ttr (Supplementary Table 2). These results
agree with our hypothesis that the upregulation of Ttr seen in TBI
is an indicator of thyroid hormone deficiency in the brain and T4
treatment reverses this change. Although our gene expression
data suggests that Ttr is more strongly modulated by T4
compared to the other known T4 transporters and receptors,
substrate binding experiments are needed to test whether Ttr is
the main T4 transporter. Through whole transcriptome profiling,
we also identified a cascade of genes and pathways potentially
involved in T4 effects. T4 treatment affected a total of 951 DEGs
involved in diverse pathways that could support cell functions
ranging from metabolic processes, cell differentiation, to cell?cell
signaling (Supplementary Data 5). Among the T4 DEGs, 121
were also altered by mTBI (Fig. 9f) and T4 reversed the
TBIinduced changes of 93 genes (Fig. 9g), including Ttr, Wdr72
(implicated in cognitive processing speed53), and Tpx2 (protective
of neurocyte apoptosis in an AD model54) (Fig. 9h). The 93 mTBI
genes normalized by T4 (Fig. 9g) were enriched for hormone
response and metabolic pathways (Fig. 9i). However, there are
pathways that are affected by TBI but are not reversed by T4
treatment as well as genes affected by T4 treatment but are not
perturbed by TBI (Supplementary Table 2; Fig. 9i), suggesting
that T4 normalizes select aspects of TBI pathology.
NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06222-0
sno 23p LAyrh6gga6pe32
High-throughput parallel single cell sequencing analysis enabled
us to characterize novel cell types, subtypes, and gene markers
within the complex cytoarchitecture of the hippocampus under
both physiological and pathological conditions, thus providing a
rich resource for the neuroscience community to study
hippocampus-related processes and diseases. In particular, our
study unveiled key aspects of the hippocampal molecular and
cellular adaptations during the acute phase of mTBI. Viewed
holistically, our analyses of cell proportion changes, global
transcriptome pattern shifts, differentially expressed genes, and
perturbed pathways indicate that the majority of hippocampal cell
types, including previously undefined cell populations,
demonstrate various degrees of sensitivity to mTBI. As the transcriptome
can instruct the functions of cells organized in circuits important
for processing of high order information, our results on the
numerous alterations in cell-type specific genes and pathways and
the reorganization of cell?cell gene co-expression patterns
provide detailed molecular information crucial for our understanding
of the mTBI pathology.
Our single cell information was also indispensable for the
prioritization of numerous potential therapeutic pathways
affected by mTBI in diverse cell types, including the thyroid hormone
pathway as indicated by Ttr. The encouraging results from the use
of T4 to normalize the disrupted pathways to improve cognitive
behavior such as learning and memory support the potential of
using single cell approaches to identify target pathways of
therapeutic applications. It is important to note that such signals
across multiple cell types can only be derived from studying all
cell types individually as done in the current study.
The transcriptome perturbations by mTBI in individual cell
types help pinpoint the cell origins of processes that likely guide
mTBI pathogenesis, such as metabolic dysfunction in astrocytes
and neurons and amyloid regulation involving ependymal and
endothelial cells. These cell-level transcriptome patterns depict
gene programs that may regulate and predict susceptibility to
post-mTBI neurological disorders such as AD, PD, PTSD,
neuroticism, and epilepsy. The information also has the potential to
guide treatments to improve mTBI outcome by targeting either
specific cell types or broad cell interactions. For instance, the fact
that genes involved in cell metabolism and amyloid processes
were recurring findings across cell types and analytical
approaches (cell-type specific DEGs, pan-hippocampal DEGs, pathways,
and cell?cell gene co-expression) suggests that these are core
processes in mTBI pathogenesis. Our single cell study opens new
avenues to deconvolute the pathogenic processes involved in
iillittffrneaoedenC irsseneuegoN euN rcee iilllllscgnaneegC iliitscygnannapgS ilittfcenanAm tIysseeunmmm
T4/TBI vs TBI (T4 unique) TBI vs sham (TBI unique)
mTBI in individual brain cell populations and prioritize both
celltype specific (such as Arhgap32 in neurons and Apoe in
ependymal cells) and potential broad-spectrum pan-hippocampal
targets (such as transthyretin and humanin).
Concussive brain injury is the most common form of brain
injury in sports, domestic, and military settings and has been
associated with numerous long-lasting and debilitating
neurological consequences. The identified genes, processes, and cell types
vulnerable to acute concussive injury will form the foundation for
mechanistic studies and for the development of novel therapeutic
strategies for mTBI and related neurological disorders. The cell
type-specific signals may also facilitate future development of
molecular diagnostic markers considering difficulties to diagnose
concussive injury using conventional neuroimaging
examinations3. Future in-depth studies to test the functionality as well as
the potential diagnostic and therapeutic values of the predicted
novel cell types, cell type specific mTBI targets, and cell-cell
coordination are warranted.
environmentally controlled rooms (22?24 ?C) with a 12 h light/dark cycle. Mice
were randomized to receive either FPI or Sham surgeries, with no investigator
blinding. FPI was performed with the aid of a microscope55 (Wild, Heerburg,
Switzerland), where a 1.5-mm diameter craniotomy was made 2.5 mm posterior to
the bregma and 2.0 mm lateral (left) of the midline with a high-speed drill (Dremel,
Racine, WI, USA). A plastic injury cap was placed over the craniotomy with
silicone adhesive and dental cement. When the dental cement hardened, the cap
was filled with 0.9% saline solution. Anesthesia was discontinued, and the injury
cap was attached to the fluid percussion device. At the first sign of hind-limb
withdrawal to a paw pinch, a mild fluid percussion pulse (1.5?1.7 atm, wake up
time greater than 5 min) was administered. Sham animals underwent an identical
preparation with the exception of the lesion. Immediately following response to a
paw pinch, anesthesia was restored and the skull was sutured. Neomycin was
applied on the suture and the mice were placed in a heated recovery chamber for
approximately an hour before being returned to their cages. After 24 h, mice were
sacrificed and fresh hippocampal tissue was dissected for use in Drop-seq (n = 3/
group with one animal per group per day; sample size was determined based on
previous single cell studies that demonstrated sufficient statistical power). All
experiments were performed in accordance with the United States National
Institutes of Health Guide for the Care and Use of Laboratory Animals and were
approved by the University of California at Los Angeles Chancellor?s Animal
Animals and mild fluid percussion injury (FPI). Male C57BL/6 J (B6) mice of
10 weeks of age (Jackson Laboratory, Bar Harbor, ME, USA) weighing between 20
and 25 g were group housed in cages (n = 3?4/group) and maintained in
Tissue dissociation for Drop-seq. The protocol by Brewer et al56. was used to
suspend cells at a final concentration of 100 cells/?l in 0.01% BSA-PBS by digesting
freshly dissected hippocampus tissue with papain (Worthington, Lakewood, NJ,
USA). Briefly, hippocampi were rapidly dissected from the ipsilateral side of the
brain on ice. The hippocampi were transferred into 4 ml HABG (Fisher Scientific,
Hampton, NH, USA) and incubated in water bath at 30 ?C for 8 min. The
supernatant was discarded and the remaining tissue was incubated with papain (12
mg in 6 ml HA-Ca) at 30 ?C for 30 min. After incubation, the papain solution was
removed from the tissue and washed with HABG three times. Using a siliconized
9in Pasteur pipette with a fire-polished tip, the solution was triturated approximately
ten times in 45 s. Next, the cell suspension was carefully applied to the top of the
prepared OptiPrep density gradient (Sigma Aldrich, St. Louis, MO, USA) and
floated on top of the gradient. The gradient was then centrifuged at 800 g for 15
min at 22 ?C. We aspirated the top 6 ml containing cellular debris. To dilute the
gradient material, we mixed the desired cell fractions with 5 ml HABG. The cell
suspension containing the desired cell fractions was centrifuged for 3 min at 22 ?C
at 200 g, and the supernatant containing the debris was discarded. Finally, the cell
pellet was loosened by flicking the tube and the cells were re-suspended in 1 ml
0.01% BSA (in PBS). This final cell suspension solution was passed through a
40micron strainer (Fisher Scientific, Hampton, NH, USA) to discard debris, followed
by cell counting.
Drop-seq single cell barcoding and library preparation. Barcoded single cells, or
STAMPs (single-cell transcriptomes attached to microparticles), and cDNA
libraries were generated following the drop seq protocol from Macosko et al12. and
version 3.1 of the online Drop-seq protocol [http://mccarrolllab.com/download/
905/]. Briefly, single cell suspensions at 100 cells/?l, EvaGreen droplet generation
oil (BIO-RAD, Hercules, CA, USA), and ChemGenes barcoded microparticles
(ChemGenes, Wilmington, MA, USA) were co-flowed through a FlowJEM
aquapel-treated Drop-seq microfluidic device (FlowJEM, Toronto, Canada) at
recommended flow speeds (oil: 15,000 ?l/hr, cells: 4000 ?l/hr, and beads 4000 ?l/hr)
to generate STAMPs. The following modifications were made to the online
published protocol to obtain enough cDNA as quantified by a high sensitivity
BioAnalyzer (Agilent, Santa Clara, CA, USA) to continue the protocol: (
) The number
of beads in a single PCR tube was 4000. (
) The number of PCR cycles was 4 + 11
) Multiple PCR tubes were pooled. The libraries were then checked on a
BioAnalyzer high sensitivity chip (Agilent, Santa Clara, CA, USA) for library
quality, average size, and concentration estimation. The samples were then
tagmented using the Nextera DNA Library Preparation kit (Illumina, San Diego, CA,
USA) and multiplex indices were added. After another round of PCR, the samples
were checked on a BioAnalyzer high sensitivity chip for library quality check before
sequencing. A cell doublet rate of 5.6% was obtained by running the microfluidic
device without the lysis buffer and counting the percentage of cell doublets through
three separate runs.
Illumina high-throughput sequencing of Drop-seq libraries. The Drop-seq
library molar concentration was quantified by Qubit Fluorometric Quantitation
(ThermoFisher, Canoga Park, CA, USA) and library fragment length was estimated
using a Bioanalyzer. Sequencing was performed on an Illumina HiSeq 2500
(Illumina, San Diego, CA, USA) instrument using the Drop-seq custom read 1B
primer (GCCTGTCCGCGGAAGCAGTGGTATCAACGCAGAGTAC) (IDT,
Coralville, IA, USA). Paired end reads were generated using custom read lengths of
24 for read 1 and 76 for read 2 and an 8 bp index read for multiplexing. Read 1
consists of the 12 bp cell barcode, followed by the 8 bp UMI, and the last 4 bp on
the read are not used. Read 2 contains the single cell transcripts.
Drop-seq data pre-processing and quality control. The fastq files of the
Dropseq sequencing data were processed to digital expression gene matrices using
Dropseq tools version 1.12 [http://mccarrolllab.com/download/922/]. The protocol
outlined in the Drop-seq alignment cookbook v1.2
followed, using default parameters. Fastq files were converted to BAM format and cell
and molecular barcodes were tagged. Reads corresponding to low quality barcodes
were removed. Next, any occurrence of the SMART adapter sequence or polyA tails
found in the reads was trimmed. These cleaned reads were converted back to fastq
format to be aligned to the mouse reference genome mm10 using STAR-2.5.0c.
After the reads were aligned, the reads which overlapped with exons were tagged
using a RefFlat annotation file of mm10. A percentage of the Chemgenes barcoded
beads which contain the UMIs and cell barcodes were anticipated to have synthesis
errors. We used the Drop-seq Tools function DetectBeadSynthesisErrors to
quantify the Chemgenes beads batch quality and estimated a bead synthesis error
rate of 5?10%, within the acceptable range. Finally, a digital gene expression matrix
for each sample was generated where each row is the read count of a gene and each
column is a unique single cell. The transcript counts of each cell were normalized
by the total number of UMIs for that cell. These values are then multiplied by
10,000 and log transformed. Digital gene expression matrices from the six samples
(3 Sham and 3 TBI samples) were combined to create three different pooled digital
gene expression matrices for: (
) all Sham samples, (
) all TBI samples and (
combined Sham and TBI samples. Single cells were identified from background
noise by using a threshold of at least 500 genes and 900 transcripts.
Identification of cell clusters. The Seurat R package version 184.108.40.206 [https://
github.com/satijalab/seurat] was used to project all sequenced cells onto two
dimensions using t-SNE and density-based spatial clustering (DBSCAN) was used
to assign clusters. To further refine the neuronal cell clusters, the BackSPIN
software14 was used to perform biclustering of the single cells identified to be neuronal
cells to further resolve this cell type. Biclustering has been previously demonstrated
to differentiate between cell types which cannot be captured by traditional
t-SNEbased approaches14. BackSPIN was run with default parameters, selecting for the
top 2000 most highly variable genes and proceeding with 5 levels of biclustering.
Identification of marker genes of individual cell clusters. We defined cell cluster
specific marker genes from our Drop-seq dataset using a bimodal likelihood ratio
test57. To determine the marker genes, the single cells were split into two groups for
each test: the cell type of interest and all remaining single cells. To be considered in
the analysis, the gene had to be expressed in at least 30% of the single cells from
one of the groups and there had to be at least a 0.25 log fold change in gene
expression between the groups. Multiple testing was corrected using the
Benjamini?Hochberg method and genes with an FDR < 0.05 are defined as marker
genes. We explored the gene-gene correlation within-group and between-group for
each cell type and confirmed the consistency of cell type identification between
samples (Supplementary Fig. 1).
Resolving cell identities of the cell clusters. We use two methods to resolve the
identities of the cell clusters. First, known cell-type specific markers from previous
studies were curated and checked for expression patterns in the cell clusters. A
cluster showing high expression levels of a known marker gene specific for a
particular cell type was considered to carry the identity of that cell type. Second, we
evaluated the overlap between known marker genes of various cell types with the
marker genes identified in our cell clusters. Overlap was assessed using a Fisher?s
exact test and significance was set to Bonferroni-corrected p < 0.05. A cluster was
considered to carry the identity of a cell type if the cluster marker genes showed
significant overlap with known markers of that cell type. The two methods showed
consistency in cell identity determination.
Known markers for major hippocampal cell types and neuronal subtypes were
retrieved from Zeisel et al.14, Habib et al.15, and the Allen Brain Atlas58. These
markers were sufficient to define all major cell types as well as GABAergic neurons,
dentate gyrus (DG), CA1 and CA3 pyramidal neurons.
Quantitative assessment of global transcriptome shifts. For each cell type, we
generate two representative cells, one for the Sham group and the other for mTBI
condition by calculating the average gene expression of each gene for each group
within that cell type. We then calculate the Euclidian distance in gene expression
between these representative cells as a metric to quantify the effect of TBI on each
cell type. To determine if the observed Euclidian distance between Sham and mTBI
cells within each cell type is significantly larger than that of random cells, we
estimated a null distribution by calculating the Euclidian distance between
randomly sampled cells of the given cell type. This permutation approach is repeated
for a total of 1000 times to generate the null distribution, which is compared to the
Euclidian distance generated from the true TBI and Sham groups to determine an
empirical p value. To correct for multiple testing across all the cell types tested, we
applied a Bonferroni correction to retrieve adjusted p values.
Cell?cell gene co-expression analysis. We assessed cell?cell gene co-expression
based on gene-level correlation patterns between any two given cell types (Fig. 6a).
To infer directionality of the interactions between two cell types, we defined a cell
type whose marker genes encode secreted peptides based on Uniprot information
as the source cell type, and then correlated the peptide-encoding marker genes
from the source cell type with genes in the target cell type. To deal with the sparse
nature of single cell data, we averaged the gene expression of Sham and TBI
samples respectively for each cell type. An interaction score is calculated from the
sum of the correlation p-values for each peptide, assuming that a peptide from a
source cell type with strong correlations with many genes in the target cell type
would have a high score and indicate strong interactions. To determine the
significance of interaction, we use a permutation approach in which a null
distribution is drawn from the interactions scores generated by the correlations between
source peptides and target genes where the expression values for each target gene
has been shuffled independently. The genes correlated with each peptide were
tested for pathway enrichment in KEGG, Reactome, BIOCARTA, GO Molecular
Functions, and GO Biological Processes to infer the key pathways involved in the
To assess the validity of our method, we benchmarked how well our gene
coexpression analysis could recapitulate the known cell?cell interactions in the
hippocampal trisynaptic circuit. As the neuronal cell types involved in this neural
circuit and present in our data are glutamatergic neurons (DG granule cells, CA1
pyramidal neurons, and CA3 pyramidal neurons) we used genes involved in
glutamate secretion to define source cell types and used the same scoring and
permutation approach above to determine significance of gene coexpression
between cell types. Due to the very small numbers of CA1 and CA3 cells in our
dataset hence limited power, we utilized data from a previous DroNc-seq study59
involving much larger populations of CA1 pyramidal cells, CA3 pyramidal cells,
and DG granule cells from 6 control animals.
Identification of DEGs between Sham and TBI. Within each identified cell type,
Sham and TBI samples are compared for differential gene expression using a
bimodal likelihood ratio test57. To be considered in the analysis, the gene had to be
expressed in at least 30% of the single cells from one of the two groups for that cell
type and there had to be at least a 0.25 log fold change in gene expression between
the groups. We correct for multiple testing using the Benjamini?Hochberg method
and genes with an FDR < 0.05 were used in downstream pathway enrichment
analyses (unless explicitly noted that a p-value of 0.01 was used instead to retrieve
suggestive pathways). Enrichment of pathways from KEGG, Reactome,
BIOCARTA, GO Molecular Functions, and GO Biological Processes was assessed with
Fisher?s exact test, followed by multiple testing correction with the
Comparison of single cell vs. bulk tissue DEGs. To define the advantages gained
by employing single cell sequencing, we simulated bulk tissue gene expression by
averaging the gene expression across all TBI single cells and all Sham single cells for
each animal. These in silico bulk tissue Sham and TBI samples were then compared
for differential gene expression using a bimodal likelihood ratio test57. FDR < 5%
was used as the cutoff to determine tissue-level DEGs, which were then compared
against those from the single-cell analysis.
RNAscope in situ hybridization. We used RNAscope Multiplex in situ
hybridization (Advanced Cell Diagnostics) to assess the expression of the genes of interest
(Wang, J Molec Diag, 2012), according to manufacturer user manual for fresh
frozen tissue. For each gene, n = 8/group was used. Two fresh frozen brain slices
from each hippocampus (15 ?m) were mounted on glass slides in 4% neutral
buffered paraformaldehyde (Fisher Scientific) for 15 min at 4 ?C. We rinsed the
slides twice in PBS, dehydrated them in 50, 70, 100, and 100% ethanol, and stored
slices in fresh 100% ethanol for overnight at ?20 ?C. Slides were dried at room
temperature for 5 min. We drew a hydrophobic barrier on slides around the brain
slices to avoid the spreading of solutions, and then treated the slides with protease
IV at room temperature for 30 min followed by a PBS wash. We applied probes to
co-localize the target gene and a cell marker gene in the same tissue to the slides
and incubated them at 40 ?C for 2 h in the HybEZ oven. Each RNAscope target
probe contains a mixture of multiple ZZ oligonucleotide probes that are bound to
the target RNA: Ttr probe (GenBank accession number NM_013697.5; target
region 141?1149, Arhgap32 probe (GenBank accession number NM_001195632.1;
target nt region, 1880?2864), Id2 probe (GenBank accession number
NM_010496.3; target nt region, 57?794), Foxj1 probe (GenBank accession number
NM_008240.3; target nt region, 1121?1836); Ly6g6e probe (GenBank accession
number NM_027366.1; target nt region, 2?902), Prox1 probe (GenBank accession
number NM_008937.2; target nt region, 590?1769), Tmem119 probe (GenBank
accession number NM_146162.2; target nt region, 2?1106), Wfs1 probe (GenBank
accession number NM_011716.2; target nt region, 757?1629), Klhl2 probe
(GenBank accession number NM_178633.3; target nt region, 534?1477), P2ry12
(GenBank accession number NM_027571.3; target nt region, 739?1854). The
negative control sections received RNase treatment before performing the
RNAscope assay, and the positive control contained a housekeeping RNA (RTU mixture
of three probes targeting POLR2A in channel C1, PPIB in channel C2, and UBC in
channel C3). We incubated the slides with amplifier probes (AMP1, 40 ?C for 30
min; AMP2, 40 ?C for 15 min; AMP3, 40 ?C for 30 min). We then incubated the
slides with fluorescently labeled probes by selecting a specific combination of colors
associated with each channel: green (Alexa 488 nm), orange (Alexa 550 nm), and
far red (Alexa 647 nm). We used AMP4 FL Alt A to detect genes of interest. We
washed the slides with 1 ? washing buffer twice in between incubations. After air
drying the slides we coverslipped them with a ProLong Gold Antifade mounting
medium containing DAPI (Invitrogen). Fluorescent images were captured at 400?
magnification (Zeiss Imager.Z1; Gottingen, DE) using the Axiovision software
(Carl Zeiss Vision, version 4.6). According to manufacturer, each dot represents a
single molecule of mRNA.
RNAscope quantification. RNAscope ISH was quantified using the automated
FISH-quant software44. Briefly, each stained tissue section was imaged in a z-stack
of 5 images with 555 nm of distance between each image. Each set of 5 z-stack
images for DAPI, Alexa 488 nm, Alexa 550 nm, and/or Alexa 647 nm depending on
the sets of genes being multiplexed were collapsed into a single.tif file for each
channel which were used in the FISH-quant software with the following settings:
pixel size?555 nm (x, y, and z), refractive index?1.0002, numeric aperture?0.75,
554, and 576 excitation and emission wavelength (Alexa 550 nm), 501 and 523
excitation and emission wavelength (Alexa 488 nm). Prior to cell segmentation, a
2D maximal projection was made from each image z-stack using a TENG focus
operator for DAPI and a HELM focus operator for Alexa 550 nm and Alexa 488
nm images with global focus measurement. These 2D projections were used in the
CellProfiler43 software to define nuclei and cell segmentations. Using these outlines
in FISH-quant, preliminary parameters for each gene were determined from a set
of 3 representative images. These parameters were applied in batch-mode to
quantify all images, and optimal threshold parameters for determining true mRNA
spots from background were defined based on distributions for: sigma-xy, sigma-z,
and the amplitude of the RNA spots. To ensure quantification of a DEG in cell
types of interest, we consider only cells which meet a certain RNA spot threshold of
the cell type specific marker gene to be the appropriate cell type. This threshold is
determined by looking for the knee of gene counts per cell. The following
thresholds were used: CA1 Pyramidal cells?6, Microglia?11, Ependymal?15, CA
Subtype 2?8, Dentate Gyrus granule cells?20. Within the group of cells which pass
this threshold, we quantify the counts per cell of the DEG of interest and compare
the count distributions between TBI and Sham groups using a bimodal likelihood
T4 treatment. L-Thyroxine sodium salt pentahydrate (T4, Sigma Chemical Co., St.
Louis, MO, USA) dissolved in saline vehicle (154 nM NaCl) was injected i.p. twice
at 1 and 6 h after FPI in the treatment group (n = 6 mice) at 1.2 ?g/100 g body
weight. Control FPI mice (n = 6) received vehicle (saline).
Behavioral tests for T4 treatment experiments. Mice from the Sham, TBI, and
T4 treatment groups were tested on the Barnes maze 7 days after injury to assess
learning acquisition and memory retention60. For learning, animals were trained
with two trials per day for 4 consecutive days, and memory retention was assessed
two days after the last learning trial. The maze was manufactured from acrylic
plastic to form a disk 1.5 cm thick and 120 cm in diameter, with 40 evenly spaced 5
cm holes at its edges. The disk was brightly illuminated (900 lumens) by four
overhead halogen lamps to provide an aversive stimulus to search for a dark escape
chamber hidden underneath a hole positioned around the perimeter of a disk. All
trials were recorded simultaneously by a video camera installed directly overhead at
the center of the maze. A trial was started by placing the animal in the center of the
maze covered under a cylindrical start chamber; after a 10 s delay, the start
chamber was raised. A training session ended after the animal had entered the
escape chamber or when a pre-determined time (5 min) had elapsed, whichever
came first. All surfaces were routinely cleaned before and after each trial to
eliminate possible olfactory cues from preceding animals. After the memory test the
animals were sacrificed immediately by decapitation and the fresh hippocampal
tissues were dissected out, frozen in liquid nitrogen, and stored in ?80 ?C for bulk
RNA-seq analysis of T4 treatment experiments. After completion of the
behavior test (11 days post-injury/treatment), hippocampal tissues were dissected
from the Sham, untreated TBI, and T4-treated TBI animals (n = 6/group for
behavior, and n = 4/group selected for sequencing). QuantSeq 3? libraries were
prepared from cDNA samples using the QuantSeq 3? mRNA-Seq Library Prep Kit
(Lexogen, Greenland, NH, USA). Libraries were run on a HiSeq4000. Reads were
aligned with STAR and read counts per gene were generated using the BlueBee
platform. Differentially expressed genes between the different groups (Sham, TBI
untreated, and TBI T4 treatment) were determined using negative binomial
models61. DEGs with p < 0.05 were included in gene signatures which were checked
for pathway enrichment. DEGs between TBI and Sham were compared for overlap
with DEGs between T4 treated mice and untreated TBI mice.
Code availability. Our analytical workflow was stored and freely accessible as a
SnakeMake file [https://github.com/darneson/DropSeq].
The sequence data that support the findings of this study have been deposited in the
Gene Expression Omnibus repository, with the series record GSE101901.
We thank Dr. Weizhe Hong for advice on neuronal subpopulation analysis and helpful
discussions. X.Y. and F.G-P. are funded by R01 DK104363 and R21 NS103088. F.GP. is
funded by R01 NS50465 and UCLA BIRC. D.A. is funded by Hyde Fellowship and
NIHNCI National Cancer Institute T32CA201160.
D.A. participated in the experimental design, collected and analyzed sequencing datasets,
and wrote the paper. Y.Z., H.R.B., I.S.A., Z.Y., G.Z conducted animal, Drop-seq, and
immunohistochemistry experiments, and edited the manuscript. X.Y. and F.G-P.
conceived the study, designed and coordinated the study, and wrote the manuscript.
Reprints and permission information is available online at http://npg.nature.com/
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.
1. Rohling , M. L. et al. A meta-analysis of neuropsychological outcome after mild traumatic brain injury: re-analyses and reconsiderations of Binder et al . ( 1997 ), Frencham et al. ( 2005 ), and Pertab et al. ( 2009 ). Clin. Neuropsychol . 25 , 608 - 623 ( 2011 ).
2. Girgis , F. , Pace , J. , Sweet , J. & Miller , J. P. Hippocampal neurophysiologic changes after mild traumatic brain injury and potential neuromodulation treatment approaches . Front. Syst. Neurosci . 10 , 8 ( 2016 ).
3. Blennow , K. et al. Traumatic brain injuries . Nat. Rev. Dis. Primers 2 , 16084 ( 2016 ).
4. Shively , S. , Scher , A. I. , Perl , D. P. & Diaz-Arrastia , R. Dementia resulting from traumatic brain injury: what is the pathology? Arch . Neurol. 69 , 1245 - 1251 ( 2012 ).
5. Lipponen , A. , Paananen , J. , Puhakka , N. & Pitk?nen , A. Analysis of posttraumatic brain injury gene expression signature reveals tubulins, Nfe2l2, Nfkb, Cd44, and S100a4 as treatment targets . Sci. Rep . 6 , 31570 ( 2016 ).
6. Meng , Q. et al. Traumatic brain injury induces genome-wide transcriptomic, methylomic, and network perturbations in brain and blood predicting neurological disorders . EBioMedicine 16 , 184 - 194 ( 2017 ).
7. Redell , J. B. et al. Analysis of functional pathways altered after mild traumatic brain injury . J. Neurotrauma 30 , 752 - 764 ( 2013 ).
8. Samal , B. B. et al. Acute response of the hippocampal transcriptome following mild traumatic brain injury after controlled cortical impact in the rat . J. Mol. Neurosci . 57 , 282 - 303 ( 2015 ).
9. von Gertten, C. , Flores Morales , A. , Holmin , S. , Mathiesen , T. & Nordqvist , A. C. Genomic responses in rat cerebral cortex after traumatic brain injury . BMC Neurosci . 6 , 69 ( 2005 ).
10. Lyeth , B. G. Historical review of the fluid-percussion TBI model . Front. Neurol . 7 , 217 ( 2016 ).
11. Mondello , S. et al. Blood-based diagnostics of traumatic brain injuries . Expert Rev. Mol. Diagn . 11 , 65 - 78 ( 2011 ).
12. Macosko , E. Z. et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets . Cell 161 , 1202 - 1214 ( 2015 ).
13. van der Maaten , L. J. P. & Hinton , G. E. Visualizing high-dimensional data using t-SNE . J. Mach. Learn. Res . 9 , 2579 - 2605 ( 2008 ).
14. Zeisel , A. et al. Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq . Science 347 , 1138 - 1142 ( 2015 ).
15. Habib , N. et al. Div-Seq: Single-nucleus RNA-Seq reveals dynamics of rare adult newborn neurons . Science 353 , 925 - 928 ( 2016 ).
16. Shlosberg , D. , Benifla , M. , Kaufer , D. & Friedman , A. Blood-brain barrier breakdown as a therapeutic target in traumatic brain injury . Nat. Rev. Neurol . 6 , 393 - 403 ( 2010 ).
17. Salehi , A. , Zhang , J. H. & Obenaus , A. Response of the cerebral vasculature following traumatic brain injury . J. Cereb. Blood Flow Metab . 37 , 2320 - 2339 ( 2017 ).
18. Lim , D. A. & Alvarez-Buylla , A. The adult ventricular-subventricular zone (VSVZ) and olfactory bulb (OB) neurogenesis . Cold Spring Harb. Perspect. Biol . 8 , a018820 ( 2016 ).
19. Carlen , M. et al. Forebrain ependymal cells are Notch-dependent and generate neuroblasts and astrocytes after stroke . Nat. Neurosci . 12 , 259 - 267 ( 2009 ).
20. Szabolcsi , V. & Celio , M. R. De novo expression of parvalbumin in ependymal cells in response to brain injury promotes ependymal remodeling and wound repair . Glia 63 , 567 - 594 ( 2015 ).
21. Butler , C. R. , Boychuk , J. A. & Smith , B. N. Effects of rapamycin treatment on neurogenesis and synaptic reorganization in the dentate gyrus after controlled cortical impact injury in mice . Front. Syst. Neurosci. 9 , 163 ( 2015 ).
22. Hayes , J. P. , Bigler , E. D. & Verfaellie , M. Traumatic brain injury as a disorder of brain connectivity . J. Int. Neuropsychol. Soc . 22 , 120 - 137 ( 2016 ).
23. Bassett , D. S. & Sporns , O. Network neuroscience . Nat. Neurosci . 20 , 353 - 364 ( 2017 ).
24. Fulcher , B. D. & Fornito , A. A transcriptional signature of hub connectivity in the mouse connectome . Proc. Natl Acad. Sci. USA 113 , 1435 - 1440 ( 2016 ).
25. Kanehisa , M. & Goto , S. KEGG: kyoto encyclopedia of genes and genomes . Nucleic Acids Res . 28 , 27 - 30 ( 2000 ).
26. Croft , D. et al. The Reactome pathway knowledgebase . Nucleic Acids Res . 42 , D472 - D477 ( 2014 ).
27. Consortium , G. O. Gene Ontology Consortium: going forward . Nucleic Acids Res . 43 , D1049 - D1056 ( 2015 ).
28. Vespa , P. et al. Metabolic crisis without brain ischemia is common after traumatic brain injury: a combined microdialysis and positron emission tomography study . J. Cereb. Blood Flow Metab . 25 , 763 - 774 ( 2005 ).
29. Van Horn , J. D. , Bhattrai , A. & Irimia , A. Multimodal imaging of neurometabolic pathology due to traumatic brain injury . Trends Neurosci . 40 , 39 - 59 ( 2017 ).
30. Sanders , M. J. , Sick , T. J. , Perez-Pinzon , M. A. , Dietrich , W. D. & Green , E. J. Chronic failure in the maintenance of long-term potentiation following fluid percussion injury in the rat . Brain Res . 861 , 69 - 76 ( 2000 ).
31. Cai , H. et al. Metabolic dysfunction in Alzheimer's disease and related neurodegenerative disorders . Curr. Alzheimer Res . 9 , 5 - 17 ( 2012 ).
32. Elliott , R. C. , Khademi , S. , Pleasure , S. J. , Parent , J. M. & Lowenstein , D. H. Differential regulation of basic helix-loop-helix mRNAs in the dentate gyrus following status epilepticus . Neuroscience 106 , 79 - 88 ( 2001 ).
33. Butovsky , O. et al. Identification of a unique TGF-?-dependent molecular and functional signature in microglia . Nat. Neurosci . 17 , 131 - 143 ( 2014 ).
34. Lu , L. , Cao , H. , Wei , X. , Li , Y. & Li , W. Iron deposition is positively related to cognitive impairment in patients with chronic mild traumatic brain injury: assessment with susceptibility weighted imaging . Biomed. Res. Int . 2015 , 470676 ( 2015 ).
35. Raz , E. et al. Brain iron quantification in mild traumatic brain injury: a magnetic field correlation study . AJNR Am. J. Neuroradiol. 32 , 1851 - 1856 ( 2011 ).
36. Nathoo , N. , Chetty , R., van Dellen, J. R. & Barnett , G. H. Genetic vulnerability following traumatic brain injury: the role of apolipoprotein E . Mol. Pathol . 56 , 132 - 136 ( 2003 ).
37. Matsuda , S. et al. The familial dementia BRI2 gene binds the Alzheimer gene amyloid-beta precursor protein and inhibits amyloid-beta production . J. Biol. Chem . 280 , 28912 - 28916 ( 2005 ).
38. Smith , D. J. et al. Genome-wide analysis of over 106 000 individuals identifies 9 neuroticism-associated loci . Mol. Psychiatry 21 , 749 - 757 ( 2016 ).
39. Soltysik-Espanola , M. et al. Characterization of Mayven, a novel actin-binding protein predominantly expressed in brain . Mol. Biol. Cell 10 , 2361 - 2375 ( 1999 ).
40. Mendez , M. F. , Owens , E. M. , Jimenez , E. E. , Peppers , D. & Licht , E. A. Changes in personality after mild traumatic brain injury from primary blast vs. blunt forces . Brain. Inj . 27 , 10 - 18 ( 2013 ).
41. Okabe , T. et al. RICS, a novel GTPase-activating protein for Cdc42 and Rac1, is involved in the beta-catenin-N-cadherin and N-methyl-D-aspartate receptor signaling . J. Biol. Chem . 278 , 9920 - 9927 ( 2003 ).
42. Ino , Y. , Gotoh , M. , Sakamoto , M. , Tsukagoshi , K. & Hirohashi , S. Dysadherin, a cancer-associated cell membrane glycoprotein, down-regulates E-cadherin and promotes metastasis . Proc. Natl Acad. Sci. USA 99 , 365 - 370 ( 2002 ).
43. Carpenter , A. E. et al. CellProfiler: image analysis software for identifying and quantifying cell phenotypes . Genome Biol . 7 , R100 ( 2006 ).
44. Mueller , F. et al. FISH-quant: automatic counting of transcripts in 3D FISH images . Nat. Methods 10 , 277 - 278 ( 2013 ).
45. Oliveira , S. M. , Ribeiro , C. A. , Cardoso , I. & Saraiva , M. J. Gender-dependent transthyretin modulation of brain amyloid-beta levels: evidence from a mouse model of Alzheimer's disease . J. Alzheimers Dis . 27 , 429 - 439 ( 2011 ).
46. Bodzioch , M. et al. Evidence for potential functionality of nuclearly-encoded humanin isoforms . Genomics 94 , 247 - 256 ( 2009 ).
47. Kim , J. et al. BRI2 (ITM2b) inhibits Abeta deposition in vivo . J. Neurosci . 28 , 6030 - 6036 ( 2008 ).
48. McLoughlin , D. M. et al. Mint2/X11-like colocalizes with the Alzheimer's disease amyloid precursor protein and is associated with neuritic plaques in Alzheimer's disease . Eur. J. Neurosci. 11 , 1988 - 1994 ( 1999 ).
49. Takamori , S. , Rhee , J. S. , Rosenmund , C. & Jahn , R. Identification of a vesicular glutamate transporter that defines a glutamatergic phenotype in neurons . Nature 407 , 189 - 194 ( 2000 ).
50. Madura , S. A. et al. Genetic variation in SLC17A7 promoter associated with response to sport-related concussions . Brain. Inj . 30 , 908 - 913 ( 2016 ).
51. Richardson , S. J. , Wijayagunaratne , R. C., D'Souza , D. G. , Darras , V. M. & Van Herck , S. L. Transport of thyroid hormones via the choroid plexus into the brain: the roles of transthyretin and thyroid hormone transmembrane transporters . Front. Neurosci. 9 , 66 ( 2015 ).
52. Li , J. et al. Thyroid hormone treatment activates protective pathways in both in vivo and in vitro models of neuronal injury . Mol. Cell. Endocrinol . 452 , 120 - 130 ( 2017 ).
53. Giddaluru , S. et al. Genetics of structural connectivity and information processing in the brain . Brain Struct. Funct . 221 , 4643 - 4661 ( 2016 ).
54. Liang , K. , Zhang , J., Yin , C. , Zhou , X. & Zhou , S. Protective effects and mechanism of TPX2 on neurocyte apoptosis of rats in Alzheimer's disease model . Exp. Ther. Med . 13 , 576 - 580 ( 2017 ).
55. Wu , A. , Molteni , R. , Ying , Z. & Gomez-Pinilla , F. A saturated-fat diet aggravates the outcome of traumatic brain injury on hippocampal plasticity and cognitive function by reducing brain-derived neurotrophic factor . Neuroscience 119 , 365 - 375 ( 2003 ).
56. Brewer , G. J. & Torricelli , J. R. Isolation and culture of adult neurons and neurospheres . Nat. Protoc . 2 , 1490 - 1498 ( 2007 ).
57. McDavid , A. et al. Modeling bi-modality improves characterization of cell cycle on gene expression in single cells . PLoS Comput. Biol . 10 , e1003696 ( 2014 ).
58. Lein , E. S. et al. Genome-wide atlas of gene expression in the adult mouse brain . Nature 445 , 168 - 176 ( 2007 ).
59. Habib , N. et al. Massively parallel single-nucleus RNA-seq with DroNc-seq . Nat. Methods 14 , 955 - 958 ( 2017 ).
60. Barnes , C. A. Memory deficits associated with senescence: a neurophysiological and behavioral study in the rat . J. Comp. Physiol. Psychol . 93 , 74 - 104 ( 1979 ).
61. Di , Y. , Schafer , D. , Cambie , J. & Chang , J. The NBP negative binomial models for assessing differential gene expression from RNA-seq . Stat. Appl. Genet. Mol. Biol . 10 , 1 - 28 ( 2011 ).
Additional information Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 018-06222-0.
Competing interests: The authors declare no competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license , visit http://creativecommons.org/ licenses/by/4.0/.
? The Author(s) 2018