Animals and mild fluid percussion injury (FPI)
Ten-week old male C57BL/6 J (B6) mice (Jackson Laboratory, Bar Harbor, ME, USA) weighing between 20-25g were housed in cages (n = 3–4/group) and maintained in 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 microscope(88) (Wild, Heerburg, Switzerland), where a 1.5-mm radius 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.4-1.6 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 or 7d, mice were sacrificed and fresh hippocampal and frontal cortex 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.
Tissue dissociation for Drop-seq
The protocol by Brewer et al(89). was used to suspend cells at a final concentration of 100 cells/μl in 0.01% BSA-PBS by digesting freshly dissected hippocampus and frontal cortex tissue with papain (Worthington, Lakewood, NJ, USA). Briefly, hippocampi from the ipsilateral side of the brain and frontal corticies were rapidly dissected on ice. The hippocampi and cortices 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 9-in 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 40-micron strainer (Fisher Scientific, Hampton, NH, USA) to discard debris, followed by cell counting. Separation of choroid plexus from hippocampus is particularly challenging due to the proximity of these two regions. Hence, the hippocampus we dissected included choroid plexus.
Peripheral blood preparation for Drop-seq
Retroorbital blood was collected in EDTA-treated collection tubes (BD Microtainer MAP Microtube, NJ, USA). Then mixed with 20* volume of ACK Lysing Buffer (Gibco, NY, USA) for 3-5 mins at room temperature to lyse red blood cells (RBCs). Leukocytes were separated by centrifugation at 300g for 5 mins at room temperature. The supernatant was removed without touching the pellet. The pellet was resuspended with cold PBS and centrifuged at 300g for 5 mins at 4 °C. The final pellet was resuspended in 0.01% BSA in PBS and filtered by a 40-micron strainer (Fisher Scientific, Hampton, NH, USA).
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 al(13). 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: (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 4000 (Illumina, San Diego, CA, USA) instrument using the Drop-seq custom read 1B primer (GCCTGTCCGCGGAAGCAGTGGTATCAACGCAGAGTAC) (IDT, Coralville, IA, USA) and PE100 reads were generated. Read 1 consists of the 12 bp cell barcode, followed by the 8 bp UMI, and the last 80 bp on the read are not used. Read 2 contains the single cell transcripts.
Drop-seq data pre-processing and quality control
Demultiplexed fastq files generated from Drop-seq were processed to digital expression gene matrices (DGEs) using Drop-seq tools version 1.13 (https://github.com/broadinstitute/Drop-seq) and dropEst(16). The workflow is available as modified version of the snakemake-based dropSeqPipe (https://github.com/Hoohm/dropSeqPipe) workflow and is available on github (https://github.com/darneson/dropSeqPipeDropEST). Briefly, fastq files were converted to BAM format and cell and molecular barcodes were tagged. Reads corresponding to low quality barcodes were removed and 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, introns, and intergenic regions were tagged using a RefFlat annotation file of mm10. To make use of reads aligning to intronic regions, which are not considered in Drop-seq tools v1.13, we used dropEst to construct digital gene expression matrices from the tagged, aligned reads where each row in the matrix is the read count of a gene and each column is a unique single cell. The count values for each cell were normalized by the total number of UMIs in that cell and then multiplied by 10,000 and log transformed. Single cells were identified from background ambient mRNA using thresholds of at least 200 genes and a maximum mitochondrial fraction of 15%.
Identification of cell clusters
The Seurat R package version 3.0.2 (https://github.com/satijalab/seurat) was used to project all sequenced cells onto two dimensions using UMAP(17) and Louvain(18) clustering was used to assign clusters. For consistent identification of cell types across different conditions (TBI vs sham) and timepoints (24-hr vs 7-day), samples were aligned using CCA(90) at the group level (timepoint + condition within a particular tissue). Specifically, the top 3000 features for each group were identified using variance stabilizing transformation and these were used to identify the top 40 CCs across the groups which were then used to find integration anchors to align the datasets. The integrated data was only used to identify and define cell types, all plots which are not explicitly designated as CCA and all downstream analyses were done on non-integrated data to retain the biological effect of timepoint and condition which is inherently removed during integration. Visualization of non-integrated data was achieved with UMAP and Louvain clustering. The optimal number of PCs used for UMAP and Louvain was determined using the Jackstraw permutation approach and a grid search of the parameters. Similarly, the density used to assign clusters was identified using a parameter grid search.
Identification of marker genes of individual cell clusters
We defined cell cluster specific marker genes from our Drop-seq dataset using the FindConservedMarkers function in Seurat across all the samples. Briefly, a Wilcoxon Rank Sum Test is run within each set of samples from a particular timepoint and condition and a meta p-value across all timepoints and conditions is computed to assess the significance of each gene as a marker for a cluster. Within each sample, the cells are split into two groups: single cells from the cell type of interest and all other single cells. To be considered in the analysis, the gene had to be expressed in at least 10% of the single cells form one of the groups and there had to be at least a 0.25 log fold change in gene expression between the groups. This process was conducted within each sample separately, and then a meta p-value was assessed from the p- values across all samples. Multiple testing was corrected using the Bonferroni method on the meta p-values and genes with an adjusted p-value < 0.05 were defined as cell type specific marker genes.
Resolving cell identities of the cell clusters
We used two methods to resolve the identities of the cell clusters. First, we used known cell-type specific markers curated from literature, single cell atlases(20, 91-93), previous studies in the hippocampus(94, 95), frontal cortex(96), and blood(97) to find distinct expression patterns in the cell clusters. A cluster showing unique expression of a known marker gene can be used to identify that cell type. We also used a classification approach leveraging the similarities between whole transcriptomes of our data and large datasets single cell datasets from the DropVIZ mouse brain atlas(20). Specifically, we obtained the single cell profiles of 113,171 hippocampus cells and 156,167 frontal cortex cells and their curated cell type labels from DropVIZ. Using the TransferData function in Seurat, we projected the PCA structure of the relevant reference datasets onto our query single cell data to classify our single cells as the most likely cell type from the annotated DropVIZ data.
Confirming absence of batch effect
To quantify the amount of batch effect present between samples within the same group (timepoint + condition) in our dataset, we leveraged the k-nearest-neighbor batch-effect test (kBET)(98). kBET interrogates the batch labels in local neighborhoods of the single cells and determines if the proportions of the batch labels in these neighborhoods differ from the global distribution. Specifically, a k-nearest neighbor matrix is constructed and 10% of the cells are selected to check the label distribution in that neighborhood. If the label distribution in the local neighborhood is similar to the global label distribution the chi-squared test does not reject the null and the batches are considered well mixed. kBET reports the average test rejection rate, however, we used the acceptance rate which is 1 – the rejection rate. In their paper, Büttner et al. noticed that kBET produced lower acceptance rates when used across an entire dataset compared to considering each cell type individually due to variations in cell-type frequencies between samples(98), thus we applied kBET to each cell type separately. To consider a cell type for testing between two samples, we required the presence of at least 15 single cells between the two samples. For each cell type we ran kBET 100 times and considered acceptance rates of >0.75 to be indicative of well-mixed batches based on the observed acceptance rates of ~0.75-0.9 for each cell type in an experiment in the kBET paper where PBMCs from eight individuals were processed in three batches and demultiplexed with demuxlet(99).
Quantitative assessment of global transcriptome shifts: Euclidean distance
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. We found that the top 1-20 highly expressed genes contributed the vast majority of the signal to this metric when considering normalized expression values. To give genes more equal weight, we transformed the expression of each gene to a z-score for each cell in the given cell type. To circumvent noise arising from lowly expressed genes, we only considered the top 1000 most highly expressed genes in 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.
Quantitative assessment of global transcriptome shifts: machine learning classifier
For each cell type, a SVM classifier was trained using the ‘caret’ library to predict sham and TBI labels using the top 1000 most highly expressed genes for that cell type, given there were at least 10 single cells per group. The model was trained 1000 times randomly sampling 70% of the data to train on with 10-fold cross validation and 3 repeats and tested on the other 30% of the data which was not seen during the training. The resulting classification accuracies at correctly predicting sham and TBI labels for held out cells across 1000 permutations were used to generate the box plots.
Quantitative assessment of global transcriptome shifts: subsample cells
For each cell type, if there were more than 100 cells per condition then they were randomly subsampled to 100 cells. DEGs were then calculated using a Wilcoxon Rank Sum Test (see below) across 1000 permutations and a meta p-value was derived using ‘minimump’ from the ‘metap’ package on the Bonferroni corrected p-values to get a stable estimate of the number of significant DEGs for each cell type.
Ligand-receptor cell-cell communication
To infer cell-cell communication, we used the CellPhoneDB(33) ligand-receptor based method. CellPhoneDB has curated 2486 interactions in the categories of protein-protein interactions, secreted and membrane proteins, and protein complexes. Based on their curated repository, CellPhoneDB predicts enriched receptor-ligand interactions between two cell types based on the expression of a receptor by one cell type and the corresponding ligand by another cell type. Only receptors and ligands which are expressed above 10% in cell type clusters are considered. To obtain a p-value for the interaction a null distribution is obtained by permuting the cluster labels of all cells and comparing the mean expression of the ligand and receptor from the cell types to the null distribution.
Identification of DEGs between Sham and TBI
Within each identified cell type, Sham and TBI samples are compared for differential gene expression using a Wilcoxon Rank Sum Test. To be considered in the analysis, the gene had to be expressed in at least 10% 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 corrected for multiple testing using Bonferroni correction and genes with an adjusted p-value < 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.
Association of DEGs with human GWAS genes of neurological disorders
Full human GWAS summary statistics of four neurological diseases (MS, epilepsy, Alzheimer’s and PD) were downloaded from the human GWAS catalog on 8-30-2019. To detect the association of cell type DEGs from TBI with human GWAS genes we used MSEA in the Mergeomics package(100). Briefly, each GWAS set was first trimmed to remove highly correlated SNPs using the Marker Dependency Filtering function with a LD50 threshold determined using the Hapmap linkage disequilibrium file for the CEU population. For each GWAS, SNPs were mapped to genes using relevant (hippocampus, frontal cortex, and blood) tissue-specific eQTL files from GTEx v8. Once mapped to genes, disease association p values for the corresponding marker were tested for enrichment in the cell type DEGs gene sets with a chi-squared-like test statistic followed by FDR estimation.
Humanin treatment
[Gly14]-Humanin (humanin, Sigma Chemical Co., St. Louis, MO, USA) dissolved in saline vehicle (154 mM NaCl) was injected i.p. twice at 1 and 6 h after FPI in the treatment group (n = 6 mice) at 40 µg/1 kg body weight. 1µg/100µl HNG was injected for the mice of average weight of 25g. Control FPI mice (n = 6) received vehicle (saline).
Behavioral tests for humanin treatment experiments
Mice from the Sham, TBI, and humanin treatment groups were trained on the Barnes maze 4 days prior to injury to facilitate learning and tested 7 days after injury to assess memory retention. For learning, animals were trained with two trials per day for 4 consecutive days, and memory retention was assessed 7 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.
Identification of genes and pathways reversed by humanin treatment
To identify genes reversed by humanin treatment, DEGs were identified for cell type between humanin treated animals and TBI animals using a Wilcoxon Rank Sum test (as described above). Genes that were significantly differentially expressed between humanin vs TBI samples and TBI vs sham samples with opposite fold change direction were considered to be reversed by the humanin treatment. Pathway enrichment analysis was conducted on these cell type DEG gene sets (as described above) to identify pathways reversed by humanin treatment.
Validation of gene expression using RNAscope
RNAscope Multiplex in situ hybridization (Advanced Cell Diagnostics, Newark, CA, USA) was conducted to evaluate the gene expression as described before(14). Briefly, the 10um brain section was mounted onto gelation-coated histological slides. The slides were fixed in pre-chilled 4% PFA for 15min at 4 ℃. The section was dehydrated in a series of ethanol followed by treatment of hydrogen peroxide for 10min and protease IV for 30min. The probes for the target gene and cell marker gene were mixed and applied to the slides with a 2-hour incubation at 40 ℃. The slides were incubated with preamplifiers, amplifiers, and dyes specific to probe channel. Finally, the sections were counterstained with DAPI and mounted with ProLong Gold Antifade mountant (Invitrogen, Carlsbad, CA, USA). The following probes were used: Mm-mt-Cytb (Cat. No. 517301); Mm-mt-Rnr1 (Cat. No. 834661); Mm-mt-Rnr2 (Cat. No. 590781); Mm-Plp1-C2 (Cat. No. 428181-C2).