Kee Wui Huang
Bernardo L. Sabatini*- Department of Neurobiology, Howard Hughes Medical Institute, Harvard Medical School, Boston, MA, United States
Viral vectors are essential tools for the study of neural circuits, with glycoprotein-deleted rabies viruses being widely used for monosynaptic retrograde tracing to map connectivity between specific cell types in the nervous system. However, the use of rabies virus is limited by the cytotoxicity and the inflammatory responses these viruses trigger. While components of the rabies virus genome contribute to its cytotoxic effects, the function of other neuronal and non-neuronal cells within the vicinity of the infected host neurons in either effecting or mitigating virally-induced tissue damage are still being elucidated. Here, we analyzed 60,212 single-cell RNA profiles to assess both global and cell-type-specific transcriptional responses in the mouse dorsal raphe nucleus (DRN) following intracranial injection of glycoprotein-deleted rabies viruses and axonal infection of dorsal raphe serotonergic neurons. Gene pathway analyses revealed a down-regulation of genes involved in metabolic processes and neurotransmission following infection. We also identified several transcriptionally diverse leukocyte populations that infiltrate the brain and are distinct from resident immune cells. Cell type-specific patterns of cytokine expression showed that antiviral responses were likely orchestrated by Type I and Type II interferon signaling from microglia and infiltrating CD4+ T cells, respectively. Additionally, we uncovered transcriptionally distinct states of microglia along an activation trajectory that may serve different functions, which range from surveillance to antigen presentation and cytokine secretion. Intercellular interactions inferred from transcriptional data suggest that CD4+ T cells facilitate microglial state transitions during the inflammatory response. Our study uncovers the heterogeneity of immune cells mediating neuroinflammatory responses and provides a critical evaluation of the compatibility between rabies-mediated connectivity mapping and single-cell transcriptional profiling. These findings provide additional insights into the distinct contributions of various cell types in mediating different facets of antiviral responses in the brain and will facilitate the design of strategies to circumvent immune responses to improve the efficacy of viral gene delivery.
Introduction
Viral vectors are widely used as tools in neuroscience for tracing, monitoring, and manipulating specific cell types and neural circuits. Neurotropic viral vectors are also being engineered for applications in clinical gene therapy, including variants of adeno-associated viruses (AAVs) that cross the blood-brain barrier with high efficiency (Chan et al., 2017). A wide variety of neurotropic viruses have been exploited as tools that have been essential for advances in basic neuroscience research, including DNA viruses such as herpes virus, AAV, and adenovirus, as well as RNA viruses such as Sindbis virus, vesicular stomatitis virus, and rabies virus (Davidson and Breakefield, 2003; Wickersham et al., 2007; Mundell et al., 2015; Ghanem and Conzelmann, 2016). Glycoprotein-deleted rabies viruses (RbVs) in particular are widely used as tools for monosynaptic retrograde tracing to label neurons that are presynaptic to a defined population of cells in the brain (Callaway and Luo, 2015). However, the utility of many viruses, including RbVs, are limited by the disruption of host cell processes and inflammatory responses that they trigger, which lead to impaired function and cytotoxicity in the circuits of interest. Understanding the various mechanisms by which cells of both the immune system and nervous system respond to these viruses may facilitate the design of improved tools and strategies to circumvent these caveats.
Neuroinflammatory processes are increasingly recognized for their importance in the etiology of neurological and psychiatric disorders. Recent studies have found significant associations between genes with known immune functions and diseases that include Alzheimer’s Disease and schizophrenia (Sekar et al., 2016; Henstridge et al., 2019; Kunkle et al., 2019; Mathys et al., 2019). Pro-inflammatory signaling and immunological perturbations in the dorsal raphe nucleus (DRN), the largest serotonergic center in the brain, have also been implicated in behavioral and mood disorders that include impulsivity, major depressive disorder, and bipolar disorder (Mahmood and Silverstone, 2001; Baumann et al., 2002; Matthews and Harrison, 2012; Howerton et al., 2014; Brisch et al., 2017). Microglia, the predominant resident immune cells of the central nervous system (CNS), express many of these disease-associated genes and undergo transcriptional changes and diversification in both acute and chronic models of inflammation, aging, and neurodegeneration (Mrdjen et al., 2018; Hammond et al., 2019; Jordão et al., 2019; Li et al., 2019). Transcriptional changes in microglia during aging are further accompanied by increased lymphocyte infiltration in the brain (Dulken et al., 2019). However, the functions of these transcriptionally distinct microglial subsets and the interactions that trigger these changes during microglial activation have yet to be elucidated.
Here, we used high-throughput single-cell RNA sequencing (scRNA-seq) and assessed changes in gene expression and cell-type composition of the mouse DRN following direct intracranial injection of glycoprotein-deleted RbVs, a procedure common to experiments using RbVs for connectivity mapping. Analysis of scRNA profiles from both RbV-injected and uninjected control mice revealed several types of infiltrating leukocytes in the DRN that are recruited by chemokines released from glial cell types. Analysis of the transcriptional changes by cell type also revealed both global and cell type-specific gene sets and pathways that underlie the antiviral defense response. Additionally, we identified transcriptionally distinct subsets of microglia, some of which minimally express canonical microglial genes, which may represent different functional states along an activation trajectory. Our study provides an in-depth assessment of the effects of RbV injection on the transcriptional state of distinct cell types in the brain. These results also provide insights into the shared and unique functions performed by various CNS cell types to mediate different facets of the immunological response.
Materials and Methods
Mice
C57BL/6J (The Jackson Laboratory, Stock #000664) were kept on a 12:12 regular light/dark cycle under standard housing conditions. All procedures were performed following protocols approved by the Harvard Standing Committee on Animal Care following guidelines described in the U.S. National Institutes of Health Guide for the Care and Use of Laboratory Animals.
Rabies Viruses
Unpseudotyped rabies viruses (B19G-SADΔG-EGFP, B19G-SADΔG-tdTomato) were generated in-house using procedures based on published protocols (Wickersham et al., 2010; Osakada and Callaway, 2013). Virions were amplified from existing stocks in three rounds of low-MOI passaging through BHK-B19G cells by transfer of filtered supernatant, with 3–4 days between passages. Cells were grown at 35°C and 5% CO2 in DMEM with GlutaMAX (Thermo Fisher Scientific, Waltham, MA, USA, #10569010) supplemented with 5% heat-inactivated FBS (Thermo Fisher Scientific, Waltham, MA, USA #10082147) and antibiotic-antimycotic (Thermo Fisher Scientific, Waltham, MA, USA #15240-062). Virions were concentrated from media from dishes containing virion-generating cells by first collecting and incubating with benzonase nuclease (1:1,000, Millipore, Kankakee, IL, USA #70664) at 37°C for 30 min, followed by filtration through a 0.22 μm PES filter. The filtered supernatant was transferred to ultracentrifuge tubes (Beckman Coulter #344058) with 2 ml of a 20% sucrose in dPBS cushion and ultracentrifuged at 20,000 RPM (Beckman Coulter SW 32 Ti rotor) at 4°C for 2 h. The supernatant was discarded and the pellet was resuspended in dPBS for 6 h on an orbital shaker at 4°C before aliquots were prepared and frozen for long-term storage at −80°C. Unpseudotyped rabies virus titers were estimated based on a serial dilution method counting infected HEK 293T cells and quantified as infectious units per ml (IU/ml).
Stereotaxic Surgeries
Mice were initially anesthetized with 5% isoflurane (80% oxygen) and maintained at 1–2.5% isoflurane after placement on the stereotaxic frame (David Kopf Instruments, Model 1900 Stereotaxic Alignment System). The scalp was cleaned and sterilized before an incision was made to expose the skull, and sterile ophthalmic ointment was applied to the eyes. For leveling the horizontal plane, a stereotaxic alignment tool (David Kopf Instruments, Model 1905) was used to zero the relative dorsoventral displacement of Bregma and Lambda, as defined in the Paxinos Brain Atlas (Paxinos and Franklin, 2001), for adjusting tilt of the anterior-posterior axis, and of two points equidistant to the left and right of Bregma for adjusting the tilt of the medial-lateral axis. Craniotomies were prepared using a mounted drill (David Kopf Instruments, Model 1911) with careful removal of the bone flap and overlying dura using forceps and a fine needle tip, and were covered with sterile 0.9% saline before and during the injection to prevent desiccation. Viruses were front-filled into a pulled glass pipette (Drummond Scientific, #5-000-2005) filled with mineral oil (Millipore Sigma, M3516) and connected to a 5 μl Hamilton syringe (Hamilton #84850) via polyethylene tubing filled with mineral oil. Glass pipettes were pulled to obtain a tip size of approximately 40–60 μm on a pipette puller (Sutter Instrument Company, P-97). Viruses were infused into target regions at approximately 100 nl/min using a syringe pump (Harvard Apparatus, #883015), and pipettes were slowly withdrawn (<10 μm/s) at least 10 min after the end of the infusion. Following wound closure, mice were placed in a cage with a heating pad until their activity was recovered before returning to their home cage. Mice were given pre- and post-operative oral carprofen (MediGel CPF, 5 mg/kg/day) as an analgesic, and monitored daily for at least 4 days post-surgery.
Stereotaxic Injection Coordinates and Volumes
All coordinates are relative to Bregma along the anterior-posterior axis and medial-lateral axis, and relative to the pial surface along the dorsoventral axis. “BL” denotes the distance between Bregma and Lambda. All injections used a straight vertical approach parallel to the DV (Z) axis. All injections were placed in the right hemisphere (positive ML values). Striatum (Str): AP = +0.40 mm, ML = ±2.45 mm, DV = −3.10 mm, 300 nl. dLGN: AP = −(2.00 * BL/4.20) mm, ML = +2.25 mm, DV = −3.00 mm, 150 nl. SN: AP = −(3.00 * BL/4.20) mm, ML = +1.32 mm, DV = −4.60 mm, 150 nl.
Histology
Mice were deeply anesthetized with isoflurane and transcardially perfused with 5–10 ml chilled 0.1 M PBS, followed by 10–15 ml chilled 4% paraformaldehyde in 0.1 M PBS. Brains were dissected out and post-fixed overnight at 4°C, followed by incubation in a storing/cryoprotectant solution of 30% sucrose and 0.05% sodium azide in 0.1 M PBS for at least 1–2 days to equilibrate. Fifty micrometer coronal slices were prepared on a freezing microtome (Leica Biosystems, SM2010 R). Fifty micrometer thick free-floating tissue sections were rinsed 3 × 5 min with 0.1 M PBS containing 0.5% Triton X-100 (PBST) before counterstaining with Neurotrace 435 (Thermo Fisher Scientific, Waltham, MA, USA N21479) at a concentration of 1:100 in 0.1 M PBS with 0.5% Triton X-100 for 1 h at room temperature. Slices were rinsed 4 × 5 min with 0.1 M PBS before they were mounted on glass slides in VectaShield mounting media (Vector Labs, H-1000). Fluorescence images were taken on an Olympus VS120 slide scanning microscope with a 10× air objective.
Single Cell Dissociation and RNA Sequencing
Identical dissociation methods, previously used and described in Huang et al. (2019), were applied to both RbV and Control groups. 8- to 10-week old C57BL/6J mice were pair-housed in a regular 12:12 light/dark cycle room before tissue collection. Mice were transcardially perfused with an ice-cold choline cutting solution containing neuronal activity blockers (110 mM choline chloride, 25 mM sodium bicarbonate, 12 mM D-glucose, 11.6 mM sodium L-ascorbate, 10 mM HEPES, 7.5 mM magnesium chloride, 3.1 mM sodium pyruvate, 2.5 mM potassium chloride, 1.25 mM sodium phosphate monobasic, 10 μM (R)-CPP, 1 μM tetrodotoxin, saturated with bubbling 95% oxygen/5% carbon dioxide, pH adjusted to 7.4 using sodium hydroxide). Brains were rapidly dissected out and sliced into 250 μm thick coronal sections on a vibratome (Leica VT1000) with a chilled cutting chamber filled with choline cutting solution. Approximately 4–5 coronal slices containing the dorsal raphe were then transferred to a chilled dissection dish containing a choline-based cutting solution for microdissection. The region containing the dorsal raphe was identified visually based on landmarks visible in unstained tissue that include the cerebral aqueduct, gray/white matter boundaries demarcating the borders of the periaqueductal gray, fiber tracts such as the medial longitudinal fasciculus and the superior cerebellar peduncle, the inferior colliculus, and the 2nd cerebellar lobule. Dissected tissue chunks were transferred to cold HBSS-based dissociation media (Thermo Fisher Scientific, Waltham, MA, USA, Cat. #14170112, supplemented to final content concentrations: 138 mM sodium chloride, 11 mM D-glucose, 10 mM HEPES, 5.33 mM potassium chloride, 4.17 mM sodium bicarbonate, 2.12 mM magnesium chloride, 0.9 mM kynurenic acid, 0.441 mM potassium phosphate monobasic, 0.338 mM sodium phosphate monobasic, 10 μM (R)-CPP, 1 μM tetrodotoxin, saturated with bubbling 95% oxygen/5% carbon dioxide, pH adjusted to 7.35 using sodium hydroxide) supplemented with an additional inhibitor cocktail (10 μM triptolide, 5 μg/ml actinomycin D, 30 μg/ml anisomycin) and kept on ice until dissections were completed. The remaining tissue was fixed in 4% paraformaldehyde in phosphate-buffered saline for histological verification. Dissected tissue chunks for each sample were pooled into a single tube for the subsequent dissociation steps. Tissue chunks were first mixed with a digestion cocktail (dissociation media, supplemented to working concentrations: 20 U/ml papain, 1 mg/ml pronase, 0.05 mg/ml DNAse I, 10 μM triptolide, 5 μg/ml actinomycin D, 30 μg/ml anisomycin) and incubated at 34°C for 90 min with gentle rocking. The digestion was quenched by adding dissociation media supplemented with 0.2% BSA and 10 mg/ml ovomucoid inhibitor (Worthington Cat. #LK003128), and samples were kept chilled for the rest of the dissociation procedure. Digested tissue was collected by brief centrifugation (5 min, 300 g), re-suspended in dissociation media supplemented with 0.2% BSA, 1 mg/ml ovomucoid inhibitor, and 0.05 mg/mL DNAse I. Tissue chunks were then mechanically triturated using fine-tip plastic micropipette tips of progressively decreasing size. The triturated cell suspension was filtered in two stages using a 70 μm cell strainer (Miltenyi Biotec Cat #130-098-462) and 40 μm pipette tip filter (Bel-Art Cat. #H136800040) and washed in two repeated centrifugation (5 min, 300 g) and re-suspension steps to remove debris before a final re-suspension in dissociation media containing 0.04% BSA and 15% OptiPrep (Sigma D1556). Cell density was calculated based on hemocytometer counts and adjusted to approximately 100,000 cells/ml. Single-cell encapsulation and RNA capture on the InDrop platform was performed at the Harvard Medical School ICCB Single Cell Core using v3 hydrogels based on previously described protocols (Zilionis et al., 2017). Suspensions were kept chilled and gently agitated until the cells flowed into the microfluidic device. Libraries were prepared and indexed following the protocols referenced above, and sequencing-ready libraries were stored at −80°C. Libraries were pooled and sequenced on an Illumina NextSeq 500 (High Output v2 kits).
Sequencing Data Processing
NGS data was processed using previously a published pipeline in Python available at: https://github.com/indrops/indrops (Klein et al., 2015). Briefly, reads were filtered by the expected structure and sorted by the corresponding library index. Valid reads were then demultiplexed and sorted by cell barcodes. Cell barcodes containing fewer than 250 total reads were discarded, and remaining reads were aligned to a reference mouse transcriptome (Ensembl GRCm38 release 87) using Bowtie 1.2.2 (m = 200, n = 1, l = 15, e = 100). For alignment, the mouse transcriptome was modified with the addition of genes from the SAD B19 rabies viruses and transgenes (B19N, B19P, B19M, B19L, EGFP, tdTomato, AmCyan1). Aligned reads were then quantified as UMI-filtered mapped read (UMIFM) counts. UMIFM counts and quantification metrics for each cell were combined into a single file sorted by the library and exported as a gunzipped TSV file.
Pre-clustering Filtering and Normalization
Analysis of the processed NGS data was performed in R version 3.4.4 using the Seurat package version 2.3.1 (Satija et al., 2015; Butler et al., 2018). A custom R script was used to combine the expression data and metadata from all libraries corresponding to a single batch, and cells with fewer than 500 UMIFM counts were removed. The expression data matrix (Genes × Cells) was filtered to retain genes with >5 UMIFM counts, and then loaded into a Seurat object along with the library metadata for downstream processing. The percentage of mitochondrial transcripts for each cell (percent.mito) was calculated and added as metadata to the Seurat object. Cells were further filtered before dimensionality reduction (Reads—min. 20,000, max. Inf; nUMI—min. 500, max. 18,000; nGene—min. 200, max. 6,000; percent.mito—min. -Inf, max. 0.1). Low-quality libraries identified as outliers on scatter plots of quality control metrics (e.g., unusually low gradient on the nGene vs. nUMI) were also removed from the dataset. Expression values were then scaled to 10,000 transcripts per cell and log-transformed. Effects of latent variables (nUMI, percent.mito, Sex) were estimated and regressed out using a GLM (ScaleData function, model.use = “linear”), and the scaled and centered residuals were used for dimensionality reduction and clustering.
Dimensionality Reduction and Batch Effect Correction
Canonical correlation analysis (CCA) was used for dimensionality reduction and mitigation of batch effects. We used 2,412 genes that were highly variable in at least two datasets to calculate canonical variates (CVs) using the RunMultiCCA function in Seurat. After inspection of the CVs, the first 21 CVs were used for subspace alignment using the AlignSubspace function to merge datasets into a single object.
Cell Clustering and Cluster Identification
Initial clustering was performed on the merged and CCA-aligned dataset using the first 21 aligned CVs. UMAP was used only for data visualization. Clustering was run using the FindClusters function using the SLM algorithm and 10 iterations. Clustering was performed at varying resolution values, and we chose a value of 2 for the resolution parameter for the initial stage of clustering. Clusters were assigned preliminary identities based on the expression of combinations of known marker genes for major cell classes and types. Low-quality cells were identified based on a combination of low gene counts, low UMIFM counts, high fraction transcripts from mitochondrial genes, and a high fraction of nuclear transcripts (e.g., Malat1, Meg3, Kcnq1ot1). These cells typically clustered together and were removed manually. Following the assignment of preliminary identities, cells were divided into data subsets as separate Seurat objects (neurons; astrocytes; ependymal cells; endothelial cells, pericytes, fibroblasts, and myocytes; immune cells; oligodendrocytes and polydendrocytes) for further sub clustering.
Subclustering
Subclustering was performed iteratively on each data subset to resolve additional cell types and subtypes. For immune cell types with proliferating cell populations (microglia, lymphocytes), cell cycle scores were calculated and regressed out using the ScaleData function in Seurat. Briefly, clustering was run at high resolution, and the resulting clusters were ordered in a cluster dendrogram built using the Ward2 method in hclust using cluster-averaged gene expression for calculating the Euclidean distance matrix. Putative doublets/multiplets were identified based on the expression of known marker genes for different cell types not in the cell subset (e.g., neuronal and glial markers). Putative doublets tended to separate from other cells and cluster together, and these clusters were removed from the dataset. Cluster separation was evaluated using the AssessNodes function and inspection of differentially expressed genes (DEGs) at each node. Clusters with poor separation, based on high OOBE scores and DE of mostly housekeeping genes, were merged to avoid over-separation of the data. The dendrogram was reconstructed after merging or removal of clusters, and the process of inspecting and merging or removing clusters was repeated until all resulting clusters could be distinguished based on a set of DEGs that we could validate separately.
Differential Expression Tests and Gene Set Enrichment Analysis (GSEA)
Tests for DE were performed using MAST version 1.4.1 (Finak et al., 2015). P-values were corrected using the Benjamini-Hochberg method and filtered a 5% false discovery rate (Q < 0.05). GSEA was performed using the fgsea package version 1.4.1 in R (Sergushichev, 2016). Genes were ordered by Z scores from MAST DE tests on either the MSigDB mouse Hallmark gene sets or Reactome pathways separately. Combined Z scores were used for most genes, and discrete component Z scores were used for genes in which the continuous component was returned as NA (e.g., gene was not expressed in one of the two comparison groups). Enrichment scores were calculated using fgsea (nperm = 100,000, maxSize = Inf). P-values were corrected in fgsea using the Benjamini-Hochberg method. Gene sets and pathways were obtained using the misgdbr package version 6.2.1.
Trajectory Inference
Trajectory inference was performed using monocle version 2.6.4 (Trapnell et al., 2014; Qiu et al., 2017). Raw count data from the Seurat microglia object was converted to a CellDataSet object using the importCDS function in monocle. Genes that were differentially expressed between microglial subclusters (Q value < 0.01, |average log2 fold change| ≥ 1) were set as the ordering genes. The minimum spanning tree was constructed using the reduceDimensions function (reduction_method = “DDRTree”, num_dim = 10, norm_method = “log”, residualModelFormula = “~BatchID + nUMI + percent.mito”, relative_expr = TRUE, scaling = TRUE).
Inference of Intercellular Interactions
Intercellular interactions were inferred using the cellphoneDB version 2.0 package in Python (Vento-Tormo et al., 2018; Efremova et al., 2019). A custom R script was used to export the single-cell gene expression data from the curated Seurat object into a counts text file and metadata text file as recommended by the developers. Only genes with human orthologs were used, and mouse gene symbols were converted to the human ortholog gene symbols before data export using data from the e!Ensembl web portal. Data was processed in cellphoneDB with statistical analysis (default parameters, iterations = 1,000, no sub-sampling). Data visualizations were made using pheatmap and ggplot2 in R based on the plotting functions provided in the cellphoneDB package.
Results
Recruitment of Circulating Leukocytes Into the DRN Following Rabies Virus Infection
Inflammatory responses were induced in the DRN by axonal infection of DRN neurons using glycoprotein-deleted RbVs of the SADΔG B19 strain. RbVs were stereotactically injected into a pair of brain regions that are both innervated by DRN serotonergic neurons, and included the striatum (Str), dorsal lateral geniculate nucleus (dLGN), nucleus accumbens (NAc), and substantia nigra (SN). Tissue containing the DRN was collected from RbV-injected animals (four mice two male, two female) 7 days post-injection (Figure 1A). Tissue chunks were dissociated into live whole-cell suspensions, and scRNA-seq libraries were prepared using the inDrop v3 platform (Klein et al., 2015; Zilionis et al., 2017). Inhibitors of neural spiking activity, transcription, and translation were included to reduce the effects of tissue dissociation on gene expression (Hrvatin et al., 2018). scRNA profiles from RbV-injected animals were analyzed together with cells collected from uninjected animals, and datasets were merged using CCA-based dataset alignment methods (Butler et al., 2018). Low-quality cells and putative multiplets were manually identified and discarded prior to analysis of differential gene expression (see “Materials and Methods” section). Our final merged dataset contained a total of 60,212 cells 20,581 cells in the RbV group (10,065 male, 10,516 female), and 39,631 cells in the Control group (17,496 male, 22,135 female). The Control dataset included single cell RNA profiles previous analyzed and described in a separate study (Huang et al., 2019). Cells were sequenced to a mean read depth of 62,061 reads/cell (min. = 20,001; median = 47,823; IQR = 43, 294; max. = 870,064), 2,603 UMIFMs (min = 501; median = 2,093; IQR = 1,849; max. = 17,997) and mean gene detection rate of 1,027 genes/cell (min. = 201; median = 896; IQR = 706; max. = 5,518). Separation of the cells by condition showed that each group was sequenced to comparable read depths with the RbV group having a higher mean (Control: 1st quartile = 30,670, median = 45,709, mean = 58,359, 3rd quartile = 69,928; RbV 1st quartile = 34,090, median = 52,542, mean = 69,191, 3rd quartile = 85,428). However, the UMIFM count (Control 1st quartile = 1,574, median = 2,338, mean = 2,856, 3rd quartile = 3,499; RbV 1st quartile = 1,082, median = 1,649, mean = 2,117, 3rd quartile = 2,605) and gene detection rates (Control 1st quartile = 741, median = 1,041, mean = 1,173, 3rd quartile = 1,446; RbV 1st quartile = 420, median = 612, mean = 745.8, 3rd quartile = 942) were lower in the RbV group compared to the Control group (Supplementary Figure S1).
Figure 1. Single-cell transcriptional analysis of responses to viral infection in the brain. (A) Experiment schematic. Unpseudotyped SADΔG B19 rabies viruses (RbV) were injected into 8–10 week old C57BL/6J mice. Each animal received a pair of injections into two different regions innervated by dorsal raphe nucleus (DRN) 5-HT neurons (nucleus accumbens (NAc) and substantia nigra (SN) in the schematic shown). Tissue containing the DRN was dissected 7 days post-injection and dissociated into whole-cell suspensions. scRNA-seq libraries were generated using the microfluidic-based inDrop platform. Age-matched uninjected animals were used as the Control group. (B) UMAP plot of merged dataset containing 60,212 cells. Control and RbV datasets were merged using canonical correlation analysis (CCA)-based dataset alignment methods. Individual points representing single cells are color-coded by cell class/type as shown in (E). (C) UMAP plot of the merged dataset with cells color-coded by experimental condition (RbV or Control). (D) Stacked bar plot showing the relative proportion of each cell type in RbV and Control groups. Cell class/type categories are color-coded following the same scheme as (B). (E) Left: dendrogram with cell class/type labels corresponding to the cluster labels in (B). Right: dot plot showing expression of example genes (columns) used to identify the major cell classes/types (rows). The color of each dot represents the average log-scaled expression of each gene across all cells in a given cluster, and the size of the dot represents the fraction of cells in the cluster in which transcripts for that gene were detected.
Cells in the merged dataset were clustered using a graph-based clustering algorithm in the CCA-aligned space (see “Materials and Methods” section). Inspection of genes enriched in each cell cluster showed that all of the resident cell types that we previously identified in the reference dataset were present in the RbV group (Figures 1B,E). However, there was a significant expansion in the proportion of immune cells in the RbV group (Figures 1C,D). This included an increase in the proportion of microglia (n = 5,949 cells, 16.9% of RbV vs. 6.2% of Control), but not resident macrophages (Res. MΦs; n = 698 cells, 1.02% RbV vs. 1.23% of Control). We also identified both myeloid and lymphoid cell clusters that were not found in the uninjected Control group. A large cluster of myeloid cells that are distinct from microglia and resident MΦs formed a significant proportion of the cells in the RbV-injected group (n = 3,575 cells, 17.2% of RbV vs. 0.11% of Control), and the proportion of lymphocytes was also greatly increased in the RbV-injected group (n = 1,714 cells, 8.25% of RbV vs. 0.04% of Control). The appearance of these non-resident immune cells is likely due to the infiltration of the brain parenchyma by circulating leukocytes, since blood was removed from the brain via transcardial perfusion before tissue collection. As an expected result of the transcardial perfusion, red blood cells were not found in the dataset. A comparison of the proportion of immune cells in groups sorted by RbV injection sites suggested that leukocyte infiltration scales with infection magnitude, which was assessed as the mean number of RbV-infected neurons labeled by each injection (Supplementary Figures S2A–E). While samples collected from mice that received a pair of RbV injections into the Str and dLGN showed an increase in the proportion of microglia compared to Control/Uninjected animals (Uninjected = 6.2%; Str-dLGN = 16.3%; NAc-SN = 17.5%), animals that received injections into SN had a larger increase in the proportion of lymphocytes (Uninjected = 0.04%; Str-dLGN = 0.75%; NAc-SN = 14.5%) and non-resident myeloid cells (Uninjected = 0.11%; Str-dLGN = 0.55%; NAc-SN = 30.9%; Supplementary Figures S2A,B). Histological analyses quantifying the number of RbV-infected neurons in the DRN and surrounding ventrolateral PAG with varying injection targets showed that injection of RbV into SN infected an order of magnitude more neurons than Str, NAc, or dLGN injections (Supplementary Figures S2C–E).
Rabies Virus Transcripts Are Detected in Both Neurons and Microglia
To identify projection neurons that are infected by RbVs, we calculated RbV gene set expression scores for each cell (see “Materials and Methods” section) and examined the distribution of RbV gene transcripts in each cell type in the RbV group (Supplementary Figure S3). As expected from the innervation of the injection sites by 5-HT neurons, the 5-HT neuron cluster (n = 275 cells in RbV group) had the highest average RbV expression score. However, the rate of detecting cells containing RbV transcripts was low. The low yield of RbV-infected neurons could be due to the following: (i) since cells were not sorted during tissue dissociation, only a small fraction of RbV-infected neurons were captured during the cell encapsulation and mRNA capture given the typical cell capture efficiency of 30–50%; (ii) RbV-infected neurons may have poor survival during the tissue digestion and dissociation and were therefore relatively depleted from the whole-cell suspensions; (iii) given the relatively low read depths and the lack of enrichment for RbV transcripts during library preparation, there are likely to be transcripts from RbV-infected cells that were not sequenced (drop-outs in scRNA-seq); and (iv) few RbV-infected neurons were labeled due to the low labeling efficiency of the viruses. Despite the low detection rate, the RbV-infected cells that were identified retained sufficient transcriptional information for clustering and assignment of cell class/type identity.
RbV transcripts were also detected in non-neuronal cells: microglia (n = 3,483 cells in RbV group), non-resident myeloid cells (n = 3,530 cells in RbV group), and astrocytes (n = 4,605 cells in RbV group) showed the next three highest average RbV expression scores and had more than 1 cell above the score threshold (red line in Supplementary Figure S3A). These cells are unlikely to have been directly infected by RbV, given the tropism of RbVs, the lack of transsynaptic spread by the replication-incompetent glycoprotein-deleted RbVs used, and the long-distance of the DRN from the injection sites. A comparison of RbV transcript counts between 5-HT neurons predicted to be the cells directly infected by RbVs, and the other three non-neuronal cell types showed that the maximum number of RbV transcript counts in these non-neuronal cells was at least an order of magnitude lower than 5-HT neurons (Supplementary Figure S3B). We hypothesize that the detection of RbV transcripts in these cells is due to their phagocytosis of mRNA-containing material released from infected neurons. We speculate that the non-zero background counts of RbV transcripts that we observed in other cells may also be due to the capture of free-floating RbV transcripts released from lysis of RbV-infected neurons in the brain and the physical disruption of some RbV-infected cells during tissue digestion and dissociation.
RNA viruses can be recognized via pattern recognition receptors that include RIG-I and RIG-I-like receptors (e.g., MDA5), which detect foreign RNA. Since RIG-I and RIG-I-like pathways are thought to be the primary means of detecting infection by RNA viruses, we calculated the per-cell expression score of genes in the KEGG RIG-I-like signaling pathway gene set (Supplementary Figure S3C). All neuron cell types had low scores for expression of genes involved in RIG-I-like signaling. Surprisingly, we did not observe a positive correlation between the RIG-I signaling gene set expression score and the RbV gene set expression score (RbV group, all cells, Pearson R = 0.02; Supplementary Figure S3D), whereas a positive correlation was observed between RbV gene counts despite the occurrence of dropouts (RbV group, all cells, Pearson R = 0.57; Supplementary Figure S3F). Genes downstream of RIG-I, such as Tmem173 (STING) were also low even in 5-HT neurons with high expression of RbV transcripts (Supplementary Figure S3E). However, RbV components, such as the P protein, which are expressed by the glycoprotein-deleted mutants that we used may inhibit interferon signaling in infected neurons (Brzózka et al., 2006; Faul et al., 2009; Scott and Nel, 2016). The low expression of interferon-stimulated genes (e.g., Isg15) and low RIG-I-like signaling gene set expression scores across neuronal types may also indicate an underlying difference in the function of neurons vs. glia in innate immunity—astrocytes and microglia, which form close associations with synapses and neurons, may serve more prominent roles in the detection of pathogens and neuronal infection due to the suppression of these pathways in neurons.
Identification of Global and Cell-Type-Specific Transcriptional Changes
To assess the population-level transcriptional responses to RbV infection, we first performed differential expression (DE) tests on simulated “bulk” RNA-seq samples separated into RbV and Control groups (Figure 2A). “Bulk” samples were simulated by averaging UMI counts for each gene across all cells in a group regardless of cell type. Many genes involved in antiviral immune responses that were expressed at low levels in the Control group were strongly up-regulated in the RbV group. These included interferon response genes such as Isg15, major histocompatibility complex (MHC) genes such as H2-Aa and H2-D1, and genes that are highly expressed by infiltrating leukocytes such as Ptprc. GSEA using the MSigDB Hallmark gene sets (Subramanian et al., 2005; Liberzon et al., 2015) showed that genes involved in type I and type II interferon responses were highly up-regulated, as well as various signaling pathways such as complement, IL-2, IL-6, and TNFα (Figure 2B). Genes involved in cell division were also up-regulated, consistent with the proliferation and clonal expansion of microglia and T cells upon activation. Our results are consistent with prior studies that have used bulk tissue profiling methods to evaluate transcriptional changes in the CNS following exposure to RbVs (Prosniak et al., 2001; Zhao et al., 2011, 2018).
Figure 2. Rabies infection induces both global and cell-type-specific transcriptional changes. (A) Scatter plot comparing averaged expression of each gene in simulated “bulk” RNA-seq for RbV and Control groups. Several of the genes most enriched in the RbV group are labeled in red. The gray line indicates the line of unity (x = y). (B) Dot plot of the top MSigDB hallmark gene sets that were significantly enriched and up-regulated in the RbV group (5% FDR, Benjamini-Hochberg correction, Normalized Enrichment Score > 0). (C) Bar plot showing the number of differentially expressed genes (DEGs) that were significantly changed (388 DEGs with Q value < 0.01 and |log2 fold-change| > 1). differential expression (DE) genes were either a part of global response programs (35 genes significant in ≥ 9 resident cell types) or cell type-specific responses (204 genes significant in 1 of 12 resident cell types). (D) The number of Reactome pathways that are either up-regulated (blue) or down-regulated (red) in the RbV group compared to controls. (E) Dot plot showing Reactome pathways that are altered in each type of resident cell. The size of each dot represents the number of genes in each cell type (columns), and the color of each dot indicates the normalized enrichment score of each pathway (rows). Pathways that were not significantly enriched (Q value ≥ 0.05, Benjamini-Hochberg correction) are not displayed.
We assessed the transcriptional changes in resident cell types by performing DE tests and GSEA separately for each cell type/class (see “Materials and Methods” section). Most DEGs that were strongly up- or down-regulated (388 genes with Q value < 0.01, absolute log2 fold-change > 1) were found to be differentially expressed in only 1 of the 12 resident cell types (Figure 2C) and were considered cell-type-specific DEGs (204 of 388 genes). Several genes were differentially expressed in 9 or more cell types and were considered part of a “global” response program (35 of 388 genes). These genes included many of the interferon-stimulated genes such as Isg15 and other genes involved in the inflammatory response. Pathway analysis using GSEA on Reactome pathway gene sets showed that the cell types varied in their responses. Neurons had the most down-regulated pathways (28 up-regulated, 148 down-regulated), whereas microglia had the highest number of up-regulated pathways (96 up-regulated, 61 down-regulated; Figure 2D). Reactome pathways that were globally up-regulated included Cytokine Signaling in Immune System and Adaptive Immune System, whereas globally down-regulated pathways were mostly related to metabolism, such as translation, oxidative phosphorylation/electron transport, lipid metabolism, and mRNA metabolism (Figure 2E). The down-regulation of these pathways may be driven by the antiviral response in an attempt to inhibit viral replication and spread, and may underlie the reduced UMIFM and gene detection rate in the RbV group relative to the Control group.
Several pathways were enriched in only a subset of cell types but not found to be cell type-specific. Genes involved in apoptosis were up-regulated in resident immune cells, fibroblasts, and astrocytes, but were down-regulated in neurons. mRNA splicing was also up-regulated in microglia and mature oligodendrocytes while being down-regulated in other glial cell types including developing oligodendrocytes. Genes involved in neurotransmission were also down-regulated in neurons and glial cell types that express neurotransmitter receptors or are involved in the regulation of synaptic transmission, including astrocytes and polydendrocytes. A reduction in neurotransmission may contribute to the host defense by limiting the spread of neurotropic viruses across synapses, and has been suggested to be induced by IFNγ signaling (Kim et al., 2002).
To identify the cell types and signaling molecules mediating the recruitment of infiltrating leukocytes, we sorted DE genes based on their gene ontology annotations to find differentially expressed cytokines and chemokines. Microglia and resident MΦs showed the highest increase in cytokine expression compared to other resident cell types (Supplementary Figure S4A). Cell types of the neurovascular unit, which include endothelial cells, astrocytes, pericytes, and fibroblasts/fibroblast-like cells, and ependymal cells at the interface with the ventricular system showed the next highest increase in cytokine production. Neurons showed the least increase in cytokine expression. Several pro-inflammatory chemokines, such as Cxcl9, Cxcl10, Ccl2, Ccl5, and Ccl7, were released by multiple cell types. Most cell types in the RbV group released a distinct set of cytokines (Supplementary Figure S4B), with fibroblasts expressing the largest set cytokines. In contrast to a recent study using a “viral déjà vu” model (Di Liberto et al., 2018), we did not detect Ccl2 expression in neurons (Supplementary Figure S4B, highlighted in red). We speculate that this may be due to a difference in the models used—our study examines changes that occur during the primary responses on the first encounter with the virus, whereas the “viral déjà vu” model investigates secondary responses.
Infiltrating Leukocytes Are Transcriptionally Diverse
To identify the types of immune cells involved in the response to RbV infection, we performed subclustering on the immune cell subset (Figure 3A, Supplementary Figure S5). Infiltrating leukocytes were separated into two main groups that were of either myeloid or lymphoid lineage. Iterative sub clustering resolved the lymphoid cell cluster into at least three distinct types. T cells (Cd3g, Cd3e) were comprised of both CD8+ effector T cells (Cd8a, Cd8b1, Tcf7, Prf1, Gzma), and a smaller number of CD4+ helper T cells (Cd4, Il2ra, Ctla4). Natural killer (NK) cells were also present in the lymphoid group and were identified by their expression of genes such as Klre1, Prf1, and Gzma. Unexpectedly, we did not detect any B cells despite their involvement in antibody production for the clearance of rabies virions (Hooper et al., 2009; Katz et al., 2017). The absence of B cells from our dataset may reflect differences in the method of inoculation (peripherally vs. direct intracranial injections), experimental time course, or the use of replication-competent rabies viruses in these other studies.
Figure 3. Type I interferon responses are induced by IFNβ produced in microglia. (A) UMAP plot of the immune cell subset. Individual points represent single cells, which are color-coded by their assigned cell-type identity. (B) UMAP feature plot showing Ifnb1 expression, in which the color of each cell represents the log-scaled UMI counts. Ifnb1 transcripts are only detected in a small subset of microglia. (C) UMAP feature plot for Ifng expression. CD4+ T cells are the primary source of Ifng. (D) Dot plot showing the expression of selected cytokine genes (columns) in each immune cell type (rows). The dot color represents the average log-scaled expression of each gene in a given cell cluster, and the size of the dot represents the fraction of cells in that cluster in which transcripts for that gene were detected.
Non-resident myeloid cells were transcriptionally heterogeneous and were comprised of several populations that were distinct from both of the resident myeloid cell types (Supplementary Figures S5A–I). The majority of non-resident myeloid cells were monocytes (Ccr2, Fn1, Plac8, Lyz2, Ly6c2lo/mid) and monocyte-derived macrophages (moMΦs; Ly6c2hi, Nr1h3, Ly6a, H2-Aa, Ms4a6d, Cd74). Several distinct clusters of dendritic cells (DCs) were also identified, which included monocyte-derived dendritic cells (moDCs; Il1b, Il1dr1, H2-Aa, Cd74, Ifitm1), conventional dendritic cells (cDCs; Ccr7, Il4i1, Cd74, Cacnb3), and plasmacytoid dendritic cells (pDCs; Ly6d, Siglech, Irf8, Runx2).
Type I Interferon Responses Are Mediated by the Release of IFNβ From Microglia
To identify the primary mediators of the antiviral transcriptional responses, we assessed the expression of interferons in each cell type. Of the genes encoding type I or II interferons, only transcripts for Ifnb1 and Ifng were detected in our dataset. Ifnb1 expression was restricted to a small subset of microglia (Figure 3B), while Ifng was expressed at high levels by CD4+ T cells, and at lower levels by CD8+ T cells (Figure 3C). Surprisingly, genes for IFNα were not detected in any cells, despite the presence of pDCs that are typically the main source of type I interferons in the periphery (Fitzgerald-Bocarsly et al., 2008). These results are consistent with previous reports of IFNβ production from microglia to limit viral spread (Drokhlyansky et al., 2017). However, our results contrast with previous reports that identified astrocytes as the primary source of IFNβ (Pfefferkorn et al., 2016), which may be due to differences in the methods used to identify IFNβ-producing cells. Given the small number of RbV-infected cells detected, we are unable to rule out the expression of IFNα/β from infected neurons.
In addition to interferons, many other cytokines were also differentially expressed between immune cell types (Figure 3D). The pro-inflammatory cytokine Il1a was expressed specifically in microglia, whereas Il1b was expressed in cDCs and moDCs. Infiltrating myeloid cells expressed high levels of several pro-inflammatory cytokines, including Ccl5 and Cxcl9. Microglia, resident MΦs, and CD4+ T cells also expressed Tnf, which may facilitate leukocyte infiltration via its effects on neurovascular cells and the blood-brain barrier (Shrestha et al., 2008; Chen et al., 2019). Other members of the TNF superfamily were also expressed in different cell types, including Tnfsf10 in resident MΦs and monocytes/monocyte-derived cells, Ltb in T cells, and Tnfsf11 in CD4+ T cells. Several anti-inflammatory cytokines were also expressed by specific cell types: Il10 was expressed specifically in CD4+ T cells, whereas the IL-1 receptor antagonist gene Il1rn was expressed by microglia and several infiltrating myeloid cell types.
Microglia Occupy Distinct States Along an Activation Trajectory
Since microglia had the highest number of up-regulated DE genes and pathways in response to the RbV infection, we performed sub clustering on the microglial subset to determine if there are transcriptionally distinct states or subtypes of activated microglia. Subclustering using algorithms based on the shared nearest-neighbors (SNN) graph divided microglia into at least seven distinct subclusters (Figure 4A). Subclusters differed in their proportion of cells from groups separated by condition (RbV Control) and injection sites (Uninjected/Str-dLGN/NAc-SN). Subclusters that had a higher proportion of cells from the Control group included subclusters I (77.0% in Uninjected; 12.5% in Str-dLGN; 10.4% in NAc-SN) and subcluster II (79.5% in Uninjected; 14.5% in Str-dLGN; 6.0% in NAc-SN). Subclusters with microglia primarily from the RbV group included subclusters V (1.0% in Uninjected; 22.4% in Str-dLGN; 76.6% in NAc-SN), VI (0.3% in Uninjected; 3.2% in Str-dLGN; 96.4% in NAc-SN), VII (0.5% in Uninjected; 12.0% in Str-dLGN; 87.5% in NAc-SN), and VIII (0.3% in Uninjected; 91.6% in Str-dLGN; 8.1% in NAc-SN). Subcluster III (13.2% in Uninjected; 69.8% in Str-dLGN; 17.0% in NAc-SN) and subcluster IV (59.3% in Uninjected; 6.5% in Str-dLGN; 34.2% in NAc-SN) showed intermediate proportions between RbV and Control groups.
Figure 4. Microglia are found in distinct transcriptional states along an activation trajectory. (A) UMAP plot showing subclusters of microglia found using shared nearest-neighbors (SNN) graph-based clustering (see “Materials and Methods” section). Points are single cells color-coded by subcluster assignment. (B–E) UMAP feature plots showing the expression of genes that are differentially expressed between microglia subclusters. Cells are color-coded by their log-scaled UMI counts for the indicated gene. (F) Dot plot showing expression of several genes (rows) that are differentially expressed between microglia subclusters [columns—color-coded according to (A)]. (G) The inferred trajectory of microglial activation. Genes that were differentially expressed between microglia subclusters were used for trajectory inference. Individual cells were arranged along a single unbranched trajectory based on their pseudo-time values. Top: cells are color-coded by subclusters, following the color scheme in (A). Cells from subclusters I-II are on the left-most end of the trajectory, while cells from subcluster VI are on the opposite end. The remaining subclusters were distributed along the middle of the inferred trajectory. Middle/Bottom: inferred graphs with cells color-coded by their expression of the indicated genes, quantified as log10(counts+0.1), which are differentially expressed between microglia sub-clusters.
DE tests between subclusters showed that they differ in expression of various genes that include those typically used as “markers” for microglia such as Cx3cr1, Csf1r, P2ry12, and Trem2 (Figures 4B–F). Subclusters I and II had the highest expression of these canonical “marker” genes, and are likely to be microglia in a “resting” or surveillance state. Expression of these canonical “marker” genes was lower in the remaining subclusters and anti-correlated with expression of interferon-stimulated genes such as Isg15 (all microglia, Cx3cr1 vs. Isg15, Pearson R = −0.56). Subclusters VII, VIII, V, and VI express progressively lower levels of Cx3cr1 (in the stated order), and instead expressed Isg15, Irf7, and Ccl2 at higher levels. Between subclusters V–VIII genes involved in different immunological processes such as antigen presentation (e.g., Cd74, H2-Aa, H2-Eb1) and cytokine secretion (e.g., Ifnb1, Tnf) are differentially expressed, suggesting that these subclusters are functionally distinct.
To assess whether subclusters V–VIII reflect distinct states or subtypes of microglia, we used trajectory inference methods based on the genes that are differentially expressed between the subclusters (Trapnell et al., 2014; Qiu et al., 2017). Branch points in the inferred trajectory/graph would suggest that these subclusters define distinct activation endpoints for microglia and diversification of activated microglia, whereas an unbranched trajectory with a single edge/path would instead suggest that the subclusters define discrete states along a continuous trajectory and describe transitional states that can be occupied through several stages of microglial activation. The graph that was constructed from the microglia transcriptomic data suggested that the subclusters lie along an unbranched trajectory (Figure 4G). Subclusters I and II that express high levels of Cx3cr1 and P2ry12 were enriched on one extreme of the inferred trajectory, and describe the “resting” or surveilling state of microglia. Subcluster VIII cells that are enriched for genes involved in antigen presentation, such as Cd74 and H2-Eb1, were densest in the middle of the trajectory where there is also an increase in expression of interferon response genes such as Isg15. Subcluster VI cells, which are enriched in the expression of cytokines such as Ccl2, Tnf, and Ifnb1, were found at the other extreme of the trajectory, and describe a secretory state. Given the higher abundance of infiltrating lymphocytes in the NAc-SN group (majority of subcluster VI microglia) compared to the Str-dLGN (majority of subcluster VIII microglia), we speculate that the transition of microglia from an antigen presentation state near the middle of the inferred activation trajectory to the secretory state may be mediated by interactions between microglia and T cells.
RbV Infection Alters the Structure of Intercellular Interactions in the DRN
To assess changes in intercellular communication between specific cell types induced by the immune response, we used cellphoneDB to predict interactions between the various cell types from our scRNA-seq data (Vento-Tormo et al., 2018; Efremova et al., 2019). Significant predicted interactions were assessed for RbV and Control groups separately (Supplementary Figure S6), and the difference between the two was used to infer changes in intercellular interactions (Figure 5). Predicted interactions among resident cell types under control conditions were highest between fibroblasts and cells of the neurovascular unit, although we anticipate that the number of predicted interactions we infer here is likely to be an underestimate since the analysis package may not include interactions mediated by neurotransmitters. A comparison of predicted interactions between RbV and Control groups suggested an overall decrease in intercellular communication between most resident cell types. The number of predicted interactions with infiltrating cell types, such as monocytes and dendritic cells, showed a large increase as expected from their relative absence in the Control group. Among the resident cell types, microglia and resident MΦs had the highest increase in the number of predicted interactions. In particular, there was an increase in the number of predicted interactions between microglia and CD4+ T cells, consistent with our hypothesis that microglia-T cell interactions may be involved in the progression of microglia along the activation trajectory. There was also an increase in the number of predicted interactions between microglia and 5-HT neurons relative to the other neuron types. We hypothesize that this may be indicative of increased signaling between microglia and RbV-infected cells.
Figure 5. The structure of intercellular interactions is altered by the immune response. Heatmap showing the change (RbV—Control) in the number of interactions between cell types resulting from RbV infection. Interactions were inferred from single-cell expression data using CellPhoneDB. Red colors indicate an increase in the number of interactions between a pair of cell types in the RbV group relative to Control, whereas blue colors indicate a decrease in the number of interactions.
Discussion
Here, we characterize both global and cell-type-specific transcriptional responses in the mouse DRN following intracranial injection of RbVs. We reveal the transcriptional diversity of both resident and infiltrating immune cells populations that mediate distinct aspects of the immune response in a region that is distal to the injection site but contains cell bodies of virally-infected neurons. Our results suggest that microglia and infiltrating T cells serve key roles in orchestrating CNS immune responses via Type I and Type II interferon signaling. We also describe transcriptionally distinct subsets of microglia that are likely to represent discrete transitional states along an activation trajectory during the progression of the immune response. Additionally, we outline the predicted changes in cell type-specific intercellular interactions that are potentially involved in different aspects of the multifaceted response, leading to the prediction that helper T cells may mediate the progression of microglia along an activation trajectory. Our study provides additional insights into the distinct immunological functions of various cell types in the brain and presents several testable models and hypotheses for experimental validation in future studies.
Transcriptional Heterogeneity Among Resident and Infiltrating Immune Cells
An advantage of using scRNA-seq in this study has been the use of transcriptome-wide RNA profiles to assign and identify cells into distinct classes, types, and states. Despite the relatively low sequencing depth that we incur with the use of high-throughput droplet-based methods, our dataset captures variations across more dimensions than “conventional” techniques such as immunolabeling and in situ hybridization. The use of transcriptome-wide RNA profiles along with careful cross-referencing of gene expression signatures in cell clusters with well-validated studies also reduces biases in cell type identification or selection that may be introduced by the choice of markers. Consistent with other recent scRNA-seq studies, our results highlight potential biases introduced when using genes such as Cx3cr1 for identifying microglia since several of these are down-regulated in the activated or disease-associated microglia (Mrdjen et al., 2018; Jordão et al., 2019; Li et al., 2019; Masuda et al., 2019).
In addition to resolving transcriptional differences between different immune cell types, our study also finds distinguishing gene expression features between distinct subclusters of microglia. Although our results suggest that these subclusters represent discrete states along a single activation trajectory, we cannot rule out the possibility of distinct activation endpoints for microglia in a branching trajectory. Possibly these differences may emerge when comparing across a variety of stimuli (e.g., “viral déjà vu,” LPS, AD transgenic models, experimental autoimmune encephalomyelitis) that may trigger different response pathways and a distinct set of activation trajectories from what we have observed in our study. Additionally, our study only examines transcriptional changes at a single time point after initial exposure to the virus, and may only capture a segment of microglial activation before the emergence of distinct endpoints. Future studies that systematically explore a full range of time points will help resolve the dynamics of these transcriptional responses. These studies will also clarify if the difference in leukocyte infiltration that we found to be correlated with the infection magnitude reflects a difference in the immune response or a change in the time course of a shared response. Deeper sequencing and datasets containing a larger collection of microglia may also help resolve finer levels of heterogeneity among lower-expressing genes.
Spatial Specificity and Inter-regional Variability in Immune Responses
While our study describes transcriptional changes in CNS resident cells, it is also plausible that immune responses may differ across brain regions. These could be due to differences in the composition of resident cells, even among closely related types. For instance, resident MΦs in the DRN, a serotonergic nucleus, express the serotonin receptor Htr1b, which was not detected in resident MΦs in the cortex, striatum, or ventral midbrain (Saunders et al., 2018; Huang et al., 2019). Interregional variation in immune signaling may also result from differences in the proximity of different locations to different neurovascular or ventricular features. The DRN is situated close to the cerebral aqueduct, as well as larger blood vessels in the ventrolateral periaqueductal gray that run along the anterior-posterior axis. Future studies comparing the responses in various regions (e.g., frontal cortex vs. DRN) may reveal shared immune mechanisms across the CNS, and may also identify specific brain regions that are particularly susceptible to immunological insults that may subsequently trigger profound and long-lasting behavioral changes. Studies profiling these responses in different brain regions from the same animal may also distinguish systemic from local or region-specific effects. Although the injection sites in our study are over 1 mm away from the DRN, previous studies have shown increases in ISG expression in the cerebellum following infection of the olfactory bulb (van den Pol et al., 2014). Since this study does not examine changes occurring at the injection sites, in regions devoid of RbV-infected neurons, or sham injected animals, we are unable to distinguish the effects of injury-induced systemic signaling from the infection-induced effects on the responses that we describe here. Nonetheless, the correlation between the magnitude of infection and abundance of infiltrating leukocytes that we observed suggests that many of the effects we observed were, at least partially, related to immune responses to viral infection of DRN neurons, rather than a response to systemic signals originating from the physical damage caused by the intracranial injection. Future studies using sham injection animals as another condition may help to dissociate the effects of injury and infection on the responses observed.
Identification of Virally Labeled Neurons With scRNA-seq
Many recent studies have used scRNA-seq to build detailed atlases of the diverse cell types that exist in the CNS. Relating the molecular profile of each cell type to its anatomical location and axonal projections remains as one of the main challenges in placing each of these cell types into neural circuits for functional studies, with the lower throughput of most anatomical tracing techniques being a major limiting factor. Methods for combining next-generation sequencing with connectivity mapping have therefore been of great interest, and several recently developed methods have used neurotropic viruses to either label cells for sorting and enrichment of transcripts from a projection-defined neuronal population (Ekstrand et al., 2014; Tasic et al., 2018), or the introduction of barcodes for reconstruction of neuronal connectivity from sequencing data (Kebschull et al., 2016; Oyibo et al., 2018). RbVs have been valuable tools in the study of neural connectivity since they can be used for cell type-specific transsynaptic retrograde tracing, which provides more specificity in the identity of the postsynaptic cell than conventional tracers (Wickersham et al., 2007). Although the results of our study caution against the use of the SAD B19 strain for functional studies (Reardon et al., 2016; Chatterjee et al., 2018), we find that cells with high expression of RbV transcripts retain sufficient transcriptional information for their classification into a specific cell type. Our results, therefore, support the feasibility of the joint use of RbVs and high-throughput sequencing methods for connectivity mapping. Although the detection of RbV-infected neurons was sparse in our dataset, limitations in the yield of RbV-labeled neurons for connectivity inference and network reconstruction can be overcome by enriching either for RbV-labeled cells, such as by FACS, or for RbV transcripts during library preparation. Future studies comparing the immune responses elicited by different viruses or virus strains used for neuroscience research will also provide crucial insights into the effects of these tools on the physiological properties of the cell types and neural circuits of interest.
Data Availability Statement
Sequencing data from rabies-injected animals generated in this study is available at NCBI (accession number: GSE136455). The control dataset from uninjected animals was described in an earlier publication (Huang et al., 2019), and is available at NCBI GEO (accession number: GSE134163).
Ethics Statement
All procedures were performed following protocols approved by the Harvard Standing Committee on Animal Care following guidelines described in the U.S. National Institutes of Health Guide for the Care and Use of Laboratory Animals.
Author Contributions
KH performed the experiments and data analysis. KH and BS designed the experiments and wrote the manuscript.
Funding
This work was supported by funding from the Howard Hughes Medical Institute (BS), National Institutes of Health (R01 MH100568 and R01 NS103226 to BS), a Harvard Brain Initiative Bipolar Disorder Seed Grant (BS), the HMS Department of Neurobiology Graduate Fellowship (KH), and the HMS Stuart H.Q. and Victoria Quan Fellowship in Neurobiology (KH).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank M. Hyun for assistance with tissue collection and scRNA-seq library preparation, N. E. Ochandarena for assistance with stereotaxic surgeries, and A.C. Philson for assistance with rabies virus production. We also thank the HMS ICCB Single Cell Core for assistance with scRNA-seq experiments on the InDrop platform; S. Hrvatin and A. Nagy (Greenberg Lab, HMS) for advice and help with scRNA-seq protocols and analysis pipelines; B. K. Lim (UCSD) and I. R. Wickersham (MIT) for advice and reagents for rabies virus production; the Bauer Core Facility at Harvard University for sequencing support; the HMS Neurobiology Imaging Facility for confocal microscopy support (P30 NS072030); J. Levasseur for assistance with animal husbandry; and L. Worth for administrative assistance; A. Saunders (HMS) for helpful discussions and feedback on our manuscript; and members of the Sabatini Lab for helpful discussions.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2020.00065/full#supplementary-material.
References
Baumann, B., Bielau, H., Krell, D., Agelink, M. W., Diekmann, S., Wurthmann, C., et al. (2002). Circumscribed numerical deficit of dorsal raphe neurons in mood disorders. Psychol. Med. 32, 93–103. doi: 10.1017/s0033291701004822
Brisch, R., Steiner, J., Mawrin, C., Krzyżanowska, M., Jankowski, Z., and Gos, T. (2017). Microglia in the dorsal raphe nucleus plays a potential role in both suicide facilitation and prevention in affective disorders. Eur. Arch. Psychiatry Clin. Neurosci. 267, 403–415. doi: 10.1007/s00406-017-0774-1
Brzózka, K., Finke, S., and Conzelmann, K.-K. (2006). Inhibition of interferon signaling by rabies virus phosphoprotein P: activation-dependent binding of STAT1 and STAT2. J. Virol. 80, 2675–2683. doi: 10.1128/JVI.80.6.2675-2683.2006
Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420. doi: 10.1038/nbt.4096
Callaway, E. M., and Luo, L. (2015). Monosynaptic circuit tracing with glycoprotein-deleted rabies viruses. J. Neurosci. 35, 8979–8985. doi: 10.1523/JNEUROSCI.0409-15.2015
Chan, K. Y., Jang, M. J., Yoo, B. B., Greenbaum, A., Ravi, N., Wu, W.-L., et al. (2017). Engineered AAVs for efficient noninvasive gene delivery to the central and peripheral nervous systems. Nat. Neurosci. 20, 1172–1179. doi: 10.1038/nn.4593
Chatterjee, S., Sullivan, H. A., MacLennan, B. J., Xu, R., Hou, Y., Lavin, T. K., et al. (2018). Nontoxic, double-deletion-mutant rabies viral vectors for retrograde targeting of projection neurons. Nat. Neurosci. 21, 638–646. doi: 10.1038/s41593-018-0091-7
Chen, A.-Q., Fang, Z., Chen, X.-L., Yang, S., Zhou, Y.-F., Mao, L., et al. (2019). Microglia-derived TNF-α mediates endothelial necroptosis aggravating blood brain-barrier disruption after ischemic stroke. Cell Death Dis. 10:487. doi: 10.1038/s41419-019-1716-9
Davidson, B. L., and Breakefield, X. O. (2003). Viral vectors for gene delivery to the nervous system. Nat. Rev. Neurosci. 4, 353–364. doi: 10.1038/nrn1104
Di Liberto, G., Pantelyushin, S., Kreutzfeldt, M., Page, N., Musardo, S., Coras, R., et al. (2018). Neurons under T cell attack coordinate phagocyte- mediated synaptic stripping. Cell 175, 458.e19–471.e19. doi: 10.1016/j.cell.2018.07.049
Drokhlyansky, E., Göz Aytürk, D., Soh, T. K., Chrenek, R., O’Loughlin, E., Madore, C., et al. (2017). The brain parenchyma has a type I interferon response that can limit virus spread. Proc. Natl. Acad. Sci. U S A 114, E95–E104. doi: 10.1073/pnas.1618157114
Dulken, B. W., Buckley, M. T., Navarro Negredo, P., Saligrama, N., Cayrol, R., Leeman, D. S., et al. (2019). Single-cell analysis reveals T cell infiltration in old neurogenic niches. Nature 571, 205–210. doi: 10.1038/s41586-019-1362-5
Efremova, M., Vento-Tormo, M., Teichmann, S. A., and Vento-Tormo, R. (2019). CellPhoneDB v2.0: inferring cell-cell communication from combined expression of multi-subunit receptor-ligand complexes. bioRxiv [Preprint].doi: 10.1101/680926
Ekstrand, M. I., Nectow, A. R., Knight, Z. A., Latcha, K. N., Pomeranz, L. E., and Friedman, J. M. (2014). Molecular profiling of neurons based on connectivity. Cell 157, 1230–1242. doi: 10.1016/j.cell.2014.03.059
Faul, E. J., Lyles, D. S., and Schnell, M. J. (2009). Interferon response and viral evasion by members of the family rhabdoviridae. Viruses 1, 832–851. doi: 10.3390/v1030832
Finak, G., McDavid, A., Yajima, M., Deng, J., Gersuk, V., Shalek, A. K., et al. (2015). MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16:278. doi: 10.1186/s13059-015-0844-5
Fitzgerald-Bocarsly, P., Dai, J., and Singh, S. (2008). Plasmacytoid dendritic cells and type I IFN: 50 years of convergent history. Cytokine Growth Factor Rev. 19, 3–19. doi: 10.1016/j.cytogfr.2007.10.006
Ghanem, A., and Conzelmann, K.-K. (2016). G gene-deficient single-round rabies viruses for neuronal circuit analysis. Virus Res. 216, 41–54. doi: 10.1016/j.virusres.2015.05.023
Hammond, T. R., Dufort, C., Dissing-Olesen, L., Giera, S., Young, A., Wysoker, A., et al. (2019). Single-cell RNA sequencing of microglia throughout the mouse lifespan and in the injured brain reveals complex cell-state changes. Immunity 50, 253.e6–271.e6. doi: 10.1016/j.immuni.2018.11.004
Henstridge, C. M., Hyman, B. T., and Spires-Jones, T. L. (2019). Beyond the neuron-cellular interactions early in Alzheimer disease pathogenesis. Nat. Rev. Neurosci. 20, 94–108. doi: 10.1038/s41583-018-0113-1
Hooper, D. C., Phares, T. W., Fabis, M. J., and Roy, A. (2009). The production of antibody by invading B cells is required for the clearance of rabies virus from the central nervous system. PLoS Negl. Trop. Dis. 3:e535. doi: 10.1371/journal.pntd.0000535
Howerton, A. R., Roland, A. V., and Bale, T. L. (2014). Dorsal raphe neuroinflammation promotes dramatic behavioral stress dysregulation. J. Neurosci. 34, 7113–7123. doi: 10.1523/JNEUROSCI.0118-14.2014
Hrvatin, S., Hochbaum, D. R., Nagy, M. A., Cicconet, M., Robertson, K., Cheadle, L., et al. (2018). Single-cell analysis of experience-dependent transcriptomic states in the mouse visual cortex. Nat. Neurosci. 21, 120–129. doi: 10.1038/s41593-017-0029-5
Huang, K. W., Ochandarena, N. E., Philson, A. C., Hyun, M., Birnbaum, J. E., Cicconet, M., et al. (2019). Molecular and anatomical organization of the dorsal raphe nucleus. Elife 8:e46464. doi: 10.7554/eLife.46464
Jordão, M. J. C., Sankowski, R., Brendecke, S. M., Sagar, Locatelli, G., Tai, Y.-H., et al. (2019). Single-cell profiling identifies myeloid cell subsets with distinct fates during neuroinflammation. Science 363:eaat7554. doi: 10.1126/science.aat7554
Katz, I. S. S., Guedes, F., Fernandes, E. R., and Dos Ramos Silva, S. (2017). Immunological aspects of rabies: a literature review. Arch. Virol. 162, 3251–3268. doi: 10.1007/s00705-017-3484-0
Kebschull, J. M., Garcia da Silva, P., Reid, A. P., Peikon, I. D., Albeanu, D. F., and Zador, A. M. (2016). High-throughput mapping of single-neuron projections by sequencing of barcoded RNA. Neuron 91, 975–987. doi: 10.1016/j.neuron.2016.07.036
Kim, I.-J., Beck, H. N., Lein, P. J., and Higgins, D. (2002). Interferon γ induces retrograde dendritic retraction and inhibits synapse formation. J. Neurosci. 22, 4530–4539. doi: 10.1523/JNEUROSCI.22-11-04530.2002
Klein, A. M., Mazutis, L., Akartuna, I., Tallapragada, N., Veres, A., Li, V., et al. (2015). Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell 161, 1187–1201. doi: 10.1016/j.cell.2015.04.044
Kunkle, B. W., Grenier-Boley, B., Sims, R., Bis, J. C., Damotte, V., Naj, A. C., et al. (2019). Genetic meta-analysis of diagnosed Alzheimer’s disease identifies new risk loci and implicates Aβ, tau, immunity and lipid processing. Nat. Genet. 51, 414–430. doi: 10.1038/s41588-019-0358-2
Li, Q., Cheng, Z., Zhou, L., Darmanis, S., Neff, N. F., Okamoto, J., et al. (2019). Developmental heterogeneity of microglia and brain myeloid cells revealed by deep single-cell RNA sequencing. Neuron 101, 207.10–223.10. doi: 10.1016/j.neuron.2018.12.006
Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004
Mahmood, T., and Silverstone, T. (2001). Serotonin and bipolar disorder. J. Affect. Disord. 66, 1–11. doi: 10.1016/s0165-0327(00)00226-3
Masuda, T., Sankowski, R., Staszewski, O., Böttcher, C., Amann, L., Scheiwe, C., et al. (2019). Spatial and temporal heterogeneity of mouse and human microglia at single-cell resolution. Nature 566, 388–392. doi: 10.1038/s41586-019-0924-x
Mathys, H., Davila-Velderrain, J., Peng, Z., Gao, F., Mohammadi, S., Young, J. Z., et al. (2019). Single-cell transcriptomic analysis of Alzheimer’s disease. Nature 570, 332–337. doi: 10.1038/s41586-019-1195-2
Matthews, P. R., and Harrison, P. J. (2012). A morphometric, immunohistochemical and in situ hybridization study of the dorsal raphe nucleus in major depression, bipolar disorder, schizophrenia, and suicide. J. Affect. Disord. 137, 125–134. doi: 10.1016/j.jad.2011.10.043
Mrdjen, D., Pavlovic, A., Hartmann, F. J., Schreiner, B., Utz, S. G., Leung, B. P., et al. (2018). High-dimensional single-cell mapping of central nervous system immune cells reveals distinct myeloid subsets in health, aging, and disease. Immunity 48, 380.e6–395.e6. doi: 10.1016/j.immuni.2018.01.011
Mundell, N. A., Beier, K. T., Pan, Y. A., Lapan, S. W., Göz Aytürk, D., Berezovskii, V. K., et al. (2015). Vesicular stomatitis virus enables gene transfer and transsynaptic tracing in a wide range of organisms. J. Comp. Neurol. 523, 1639–1663. doi: 10.1002/cne.23761
Osakada, F., and Callaway, E. M. (2013). Design and generation of recombinant rabies virus vectors. Nat. Protoc. 8, 1583–1601. doi: 10.1038/nprot.2013.094
Oyibo, H., Cao, C., Ferrante, D. D., Zhan, H., Koulakov, A. A., Enquist, L., et al. (2018). A computational framework for converting high-throughput DNA sequencing data into neural circuit connectivity. bioRxiv [Preprint]. doi: 10.1101/244079
Paxinos, G., and Franklin, K. B. J. (2001). The Mouse Brain in Stereotaxic Coordinates. San Diego, CA: Academic Press.
Pfefferkorn, C., Kallfass, C., Lienenklaus, S., Spanier, J., Kalinke, U., Rieder, M., et al. (2016). Abortively infected astrocytes appear to represent the main source of interferon β in the virus-infected brain. J. Virol. 90, 2031–2038. doi: 10.1128/jvi.02979-15
Prosniak, M., Hooper, D. C., Dietzschold, B., and Koprowski, H. (2001). Effect of rabies virus infection on gene expression in mouse brain. Proc. Natl. Acad. Sci. U S A 98, 2758–2763. doi: 10.1073/pnas.051630298
Qiu, X., Mao, Q., Tang, Y., Wang, L., Chawla, R., Pliner, H. A., et al. (2017). Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14, 979–982. doi: 10.1038/nmeth.4402
Reardon, T. R., Murray, A. J., Turi, G. F., Wirblich, C., Croce, K. R., Schnell, M. J., et al. (2016). Rabies virus CVS-N2c(ΔG) strain enhances retrograde synaptic transfer and neuronal viability. Neuron 89, 711–724. doi: 10.1016/j.neuron.2016.01.004
Satija, R., Farrell, J. A., Gennert, D., Schier, A. F., and Regev, A. (2015). Spatial reconstruction of single-cell gene expression data. Nat. Biotechnol. 33, 495–502. doi: 10.1038/nbt.3192
Saunders, A., Macosko, E. Z., Wysoker, A., Goldman, M., Krienen, F. M., de Rivera, H., et al. (2018). Molecular diversity and specializations among the cells of the adult mouse brain. Cell 174, 1015.e16–1030.e16. doi: 10.1016/j.cell.2018.07.028
Scott, T. P., and Nel, L. H. (2016). Subversion of the immune response by rabies virus. Viruses 8:E231. doi: 10.3390/v8080231
Sekar, A., Bialas, A. R., de Rivera, H., Davis, A., Hammond, T. R., Kamitaki, N., et al. (2016). Schizophrenia risk from complex variation of complement component 4. Nature 530, 177–183. doi: 10.1038/nature16549
Sergushichev, A. A. (2016). An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. bioRxiv [Preprint]. doi: 10.1101/060012
Shrestha, B., Zhang, B., Purtha, W. E., Klein, R. S., and Diamond, M. S. (2008). Tumor necrosis factor α protects against lethal West Nile virus infection by promoting trafficking of mononuclear leukocytes into the central nervous system. J. Virol. 82, 8956–8964. doi: 10.1128/jvi.01118-08
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U S A 102, 15545–15550. doi: 10.1073/pnas.0506580102
Tasic, B., Yao, Z., Graybuck, L. T., Smith, K. A., Nguyen, T. N., Bertagnolli, D., et al. (2018). Shared and distinct transcriptomic cell types across neocortical areas. Nature 563, 72–78. doi: 10.1038/s41586-018-0654-5
Trapnell, C., Cacchiarelli, D., Grimsby, J., Pokharel, P., Li, S., Morse, M., et al. (2014). The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386. doi: 10.1038/nbt.2859
van den Pol, A. N., Ding, S., and Robek, M. D. (2014). Long-distance interferon signaling within the brain blocks virus spread. J. Virol. 88, 3695–3704. doi: 10.1128/jvi.03509-13
Vento-Tormo, R., Efremova, M., Botting, R. A., Turco, M. Y., Vento-Tormo, M., Meyer, K. B., et al. (2018). Single-cell reconstruction of the early maternal-fetal interface in humans. Nature 563, 347–353. doi: 10.1038/s41586-018-0698-6
Wickersham, I. R., Lyon, D. C., Barnard, R. J. O., Mori, T., Finke, S., Conzelmann, K.-K., et al. (2007). Monosynaptic restriction of transsynaptic tracing from single, genetically targeted neurons. Neuron 53, 639–647. doi: 10.1016/j.neuron.2007.01.033
Wickersham, I. R., Sullivan, H. A., and Seung, H. S. (2010). Production of glycoprotein-deleted rabies viruses for monosynaptic tracing and high-level gene expression in neurons. Nat. Protoc. 5, 595–606. doi: 10.1038/nprot.2009.248
Zhao, P., Liu, S., Zhong, Z., Jiang, T., Weng, R., Xie, M., et al. (2018). Analysis of expression profiles of long noncoding RNAs and mRNAs in brains of mice infected by rabies virus by RNA sequencing. Sci. Rep. 8:11858. doi: 10.1038/s41598-018-30359-z
Zhao, P., Zhao, L., Zhang, T., Qi, Y., Wang, T., Liu, K., et al. (2011). Innate immune response gene expression profiles in central nervous system of mice infected with rabies virus. Comp. Immunol. Microbiol. Infect. Dis. 34, 503–512. doi: 10.1016/j.cimid.2011.09.003
Keywords: scRNA-seq, neuroinflammation, G-deleted rabies virus, microglia, brain infiltration
Citation: Huang KW and Sabatini BL (2020) Single-Cell Analysis of Neuroinflammatory Responses Following Intracranial Injection of G-Deleted Rabies Viruses. Front. Cell. Neurosci. 14:65. doi: 10.3389/fncel.2020.00065
Received: 16 January 2020; Accepted: 04 March 2020;
Published: 20 March 2020.
Edited by:
Albert Giralt, University of Barcelona, SpainReviewed by:
Daniel Tornero, University of Barcelona, SpainJose P. Lopez-Atalaya, Spanish National Research Council, Spain
Copyright © 2020 Huang and Sabatini. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Bernardo L. Sabatini, YnNhYmF0aW5pQGhtcy5oYXJ2YXJkLmVkdQ==