Single cell molecular alterations reveal target cells and pathways of concussive brain injury

Nature Communications, Sep 2018

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

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 Xia Yang 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. - ;:, )( 0 9 8 7 6 5 4 3 2 1 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 neuroscience. 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 discovery11. 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. Results 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 Fig. 1). 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 a 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 these results. 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 dysfunction21. 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. Pex5l Cpne4 Fam163b Slc17a6 Rarres2 Slc13a4 Arhgap12 Ly6e Stxbp6 Nts Odf3b Mgp Fibcd1 Npy2r Igfbp5 Cbln2 Dnali1 Pcolce TBI DG granule cells Neuronal subtype 2 CA1 Neurons Sham DG granule cells CA3 neurons GABAergic interneurons Max Min c d e CA1 neurons: Wfs1 DG granule cells: Dsp CA subtype 2 neurons GABAergic interneurons: Gad2 CA subtype 1 neurons Neuronal subtype 1 CA1 neurons (139) 15 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 ( 48 ) 0 2 0 0 1 0 0 0 1 0 1 2 4 6 1 0 0 0 2 2 0 0 1 1 0 0 9 5 7 4 1 1 1 1 0 12 6 3 5 5 4 2 1 0 9 6 8 4 2 3 1 1 0 0 0 0 0 0 1 0 14 11 5 1 2 3 41 35 2 3 3 0 0 0 5 0 1 31 5 4 ) CA1ABA(49 ) CA2ABA(53 ) CA3ABA(52 ) CA1divSeq(101 ) CA2divSeq(153 ) CA3divSeq(139 ) DGABA(50 ) divSeq(164 DG ) GABAergicdivSeq(122 Atp2b1,Hpca,Spink8,Arpc2,Ptn Ly6e,Atp1b1,Slc24a2,Prkca,Cpe Dcn,Nnat,Timp2,Cpne7,Ndufc2 Serpini1,Ly6g6e,Cbln2,Fezf2,Nxph3 Cplx2,Ncdn,Olfm1,Synpr,Sem5a C1ql2,Tpt1,Plekha2,Prkce,Nrn1 Gad2,Gad1,Slc32a1,Nap1l5,Slc6a1 Tshz2,Arpp21,Ddit4l,Nrep,Meis2 Ndnf,Diablo,Cacna2d2,Marcks,1500016L03Rik Corrected P value (0.1,1] (1e?18,1e?12] (0.001,0.1] (1e?25,1e?18] (1e?06,0.001] (1e?33,1e?25] (1e?09,1e?06] (1e?42,1e?33] (1e?12,1e?09] (1e?52,1e?42] f CA1 neurons CA3 neurons CA subtype 1 neurons CA subtype 2 neurons DG granule cells GABAergic interneurons 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 mTBI pathogenesis. 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 potentiation post-TBI30. 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 p U n w o D f Whole tissue Single cell 280 272 345 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 cell types. 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 p U n w o D Max Min a *** ** a 2 m t I 5 d y x F 7 k A e o p A 2 1 y r 2 P f r T 2 l h l K 2 3 p a g h r A 2 d I *** ** *** l a il e h t o d n E s n o r u e n 1 A C ** s n o r u e n 2 e p y t b u S A C *** s ll e c e l u n a r G G D *** s n o r u e n 1 A C ** * TBI *** * * * * * s ll e c e l u n a r g G D 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 DAPI P2ry12 Tmem119 DAPI sno 23p LAyrh6gga6pe32 r a u g ne rh 2TSAC -)(1AAC Discussion 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. 3 2 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 Condition Sham/veh TBI/veh TBI/T4 a 80 ) (s 60 y cn40 e t La20 0 e h ge 0.4 n a h cd 0.0 l o f g2?0.4 o L ge 4 n a h c ld 0 o f 2 og?4 L ?8 *** * TBIvsSham T4vsTBI *** Ptgds TBIvsSham T4vsTBI *** *** *** ** * ** Ttr Ucp2 Cdk8 Tpx2 * * ** * Wdr72 * * * * f TBIvsSham T4vsTBI 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 Research Committee. Methods 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 [ 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: ( 1 ) The number of beads in a single PCR tube was 4000. ( 2 ) The number of PCR cycles was 4 + 11 cycles. ( 3 ) 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 []. The protocol outlined in the Drop-seq alignment cookbook v1.2 [] was 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: ( 1 ) all Sham samples, ( 2 ) all TBI samples and ( 3 ) 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 [https://] 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 interactions. 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 Benjamini?Hochberg method. 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 ratio test57. 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 tissue RNA-sequencing. 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 []. Data availability The sequence data that support the findings of this study have been deposited in the Gene Expression Omnibus repository, with the series record GSE101901. Acknowledgements 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. Author contributions 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 reprintsandpermissions/ 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 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 licenses/by/4.0/. ? The Author(s) 2018

This is a preview of a remote PDF:

Douglas Arneson, Guanglin Zhang, Zhe Ying, Yumei Zhuang, Hyae Ran Byun, In Sook Ahn, Fernando Gomez-Pinilla, Xia Yang. Single cell molecular alterations reveal target cells and pathways of concussive brain injury, Nature Communications, 2018, DOI: 10.1038/s41467-018-06222-0