key: cord-0953485-rre4wuia
authors: Zhang, Yachun; Xing, Xudong; Long, Ben; Cao, Yandi; Hu, Simeng; Li, Xiangning; Yu, Yalan; Tian, Dayong; Sui, Baokun; Luo, Zhaochen; Liu, Wei; Lv, Lei; Wu, Qiong; Dai, Jinxia; Zhou, Ming; Han, Heyou; Fu, Zhen F.; Gong, Hui; Bai, Fan; Zhao, Ling
title: A spatial and cellular distribution of rabies virus infection in the mouse brain revealed by fMOST and single‐cell RNA sequencing
date: 2022-01-20
journal: Clin Transl Med
DOI: 10.1002/ctm2.700
sha: dfc6b70df870f2bf2850f58ee7ad131eea6d4660
doc_id: 953485
cord_uid: rre4wuia

BACKGROUND: Neurotropic virus infection can cause serious damage to the central nervous system (CNS) in both humans and animals. The complexity of the CNS poses unique challenges to investigate the infection of these viruses in the brain using traditional techniques. METHODS: In this study, we explore the use of fluorescence micro‐optical sectioning tomography (fMOST) and single‐cell RNA sequencing (scRNA‐seq) to map the spatial and cellular distribution of a representative neurotropic virus, rabies virus (RABV), in the whole brain. Mice were inoculated with a lethal dose of a recombinant RABV encoding enhanced green fluorescent protein (EGFP) under different infection routes, and a three‐dimensional (3D) view of RABV distribution in the whole mouse brain was obtained using fMOST. Meanwhile, we pinpointed the cellular distribution of RABV by utilizing scRNA‐seq. RESULTS: Our fMOST data provided the 3D view of a neurotropic virus in the whole mouse brain, which indicated that the spatial distribution of RABV in the brain was influenced by the infection route. Interestingly, we provided evidence that RABV could infect multiple nuclei related to fear independent of different infection routes. More surprisingly, our scRNA‐seq data revealed that besides neurons RABV could infect macrophages and the infiltrating macrophages played at least three different antiviral roles during RABV infection. CONCLUSION: This study draws a comprehensively spatial and cellular map of typical neurotropic virus infection in the mouse brain, providing a novel and insightful strategy to investigate the pathogenesis of RABV and other neurotropic viruses.

Our fMOST data provided the 3D view of a neurotropic virus in the whole mouse brain, which indicated that the spatial distribution of RABV in the brain was influenced by the infection route. Interestingly, we provided evidence that RABV could infect multiple nuclei related to fear independent of different infection routes. More surprisingly, our scRNA-seq data revealed that besides neurons RABV could infect macrophages and the infiltrating macrophages played at least three different antiviral roles during RABV infection.

This study draws a comprehensively spatial and cellular map of typical neurotropic virus infection in the mouse brain, providing a novel and insightful strategy to investigate the pathogenesis of RABV and other neurotropic viruses.

fear, fMOST technology, macrophages, NK cells, rabies virus, single-cell RNA-seq, whole brain distribution

Neurotropic virus infection can cause an array of immediate and delayed neuropathology in humans and animals, and nearly half of the emerging viruses can invade the central nervous system (CNS). 1 These infections pose a major challenge to humans and animal healthcare worldwide. The complex structures and functions of the CNS and the diversity of neurotropic viruses make it difficult to find an effective treatment for these diseases. By combining human pathological data with experimental animal models, virologists have much advanced our understanding of the mechanisms underlying how viruses enter the CNS and cause neurological disease, but more in-depth studies are still in urgent need to facilitate the development of novel vaccines and antiviral therapeutics. 2 As a typical neurotropic virus, rabies virus (RABV) belongs to the family Rhabdoviridae. RABV hijacks the cellular transport machinery and moves along microtubules to the nearest sensory neurons by retrograde axonal transport. 3 After replication, RABV travels along the corticospinal tract and invades into the brain where the virus can efficiently replicate and infect in most regions, resulting in fatal encephalitis. 3 However, due to the complexity of the CNS, the exact regions and cell types infected by RABV are still kept elusive, impeding the endeavour of researchers to find the therapeutic target for RABV. For many decades, light microscopy was the key tool for investigating the invasion and spread of neurotropic viruses in the host's brain. However, observing a whole animal brain, which is centimetre-sized, has been faced a lot of difficulty. Histological sections were usually prepared for observing the internal microstructures of large specimens, but it is difficult to align the serial figures of continuous ultrathin sectioning on large specimens by a traditional light microscopy. In 2010, a micro-optical sectioning tomography (MOST) system was developed, which allowed the mapping of a whole mouse brain to the single neuron level. 4 In 2013, a fluorescence MOST (fMOST) was developed by using a resin-embedding strategy to maintain fluorescence and an automated fluorescence MOST system to continually imaging. 5 Recently, a modified fMOST, benefiting from simple sample preparation and high-throughput imaging, was applied to map neural circuits. 6, 7 This fMOST system can efficiently provide wide-field optical-sectioning imaging of samples embedded in agarose. Single-cell sequencing is an emerging technology that allows high-throughput sequencing analyses of genome, transcriptome, and methylome at the singlecell level. It has been recently implemented in the study of influenza virus, 8 HIV, 9 Ebola virus 10 and SARS-CoV-2 11 infections, providing an unbiased and comprehensive visualization of the target cells that virus can infect and a global immunological responses of the host.

In this study, we constructed a recombinant RABV expressing EGFP (RABV-EGFP) and utilized the fMOST system to portray the spatial distribution of RABV in the brain of infected mice. Meanwhile, scRNA-seq was applied to pinpoint the cellular distribution of RABV infection. These two combined techniques reveal a three-dimension (3D) distribution of RABV in the whole brain and identify several major cell types infected by RABV, especially macrophages and NK cells. Our results provide a better understanding of the pathogenesis of RABV s and shed new light on the development of novel therapies for rabies. 

Groups of female C57BL/6 mice in 6-week-old were challenged with 10×LD 50 RABV-EGFP by intramuscular (i.m.), otic subcutaneous (o.s.), or intranasal (i.n.) route (n = 3). In severe paralysis, mice were anaesthetized with ketamine/xylazine and then intracardially perfused by PBS followed by 10% neutral-buffered formalin. Brains were removed and fixed with 10% neutral-buffered formalin for 24 h at 4 • C. Then each brain was rinsed with PBS overnight and embedded with oxidized agarose before imaging.

The agarose-embedded brains were imaged by fMOST system. 7 During the process of imaging, the mouse brain was immersed in water and observed with a water immersion objective (1.0 NA, 20×). The fluorescent signals were captured by a scientific complementary metaloxide semi-conductor (sCMOS) camera with high sensitivity and speed. A 3D stage moved the mouse brain mosaic-bymosaic to extend the viewing field in the imaging part, and then the stage carried the brain towards the oscillatory blade to remove the imaged tissue in the vibrating sectioning part. The full-volumetric imaging was performed with the cycle of imaging and sectioning until the wholebrain dataset was collected. For a single mouse brain, the dataset including approximately 3000 coronal sections was collected in three days at a voxel resolution of 0.32 μm × 0.32 μm ×5 μm.

Image preprocessing for mosaic stitching and uneven lateral illumination correction was performed as reported previously. 6 For whole-brain visualization of RABV infection in different routes, the RABV-infected cells were automatically identified by NeuroGPS and registered into the Allen CCFv3 as previously described. 12, 13 Briefly, we manually segmented the brain outline and several anatomical regions as landmarks, and then performed affine transformation and symmetric image normalization to acquire transformation parameters. Next, we warped the RABVinfected cells' coordinates to Allen CCFv3 using the transformation parameters for registration. The 3D visualization of brain regions and RABV distribution was visualized by Amira software (v6.1.1, FEI, France) to generate the figures of maximum intensity projection, volume and surface rendering.

To quantify the RABV-infected cells, the EGFP-labelled cells were automatically segmented and coregistered to Mouse Reference Atlas as previously described. 14, 15 Briefly, each coronal section of the brain was background subtracted, Gaussian filtered, and threshold segmented to binary image, and RABV-infected cells were segmented with individually adjusted binary thresholds according to varying fluorescent intensities. And then the soma coordinates were warped and coregistered to the corresponding Allen Mouse Reference Atlas (Allen Institute for Brain Science) coordinate using non-rigid registration with freeform deformation. The anatomical sub-regions to which infected cells belonged were then mapped to the Allen Mouse Reference Atlas. We calculated the total number and the density of EGFP-labelled infected cells for each anatomical sub-region with different infection routes (n = 3 per infection route).

Groups of C57BL/6 mice in 6-week-old were i.m. infected with 10×LD 50 RABV-EGFP in the hind limb. At the stage of paralysis and moribund, mice were anaesthetized with ketamine/xylazine and then brains were collected. Singlecell suspensions were obtained with adult brain dissociation kit (Miltenyi Biotec, catalogue # 130-107-677). EGFP positive cells were sorted using the GFP channels of the Bio-Rad S3e instrument. cDNA libraries were prepared from single-cell suspensions following the instruction of 10×Genomics 3′ V3. Cells obtained from the whole brain of the healthy mice with the same method and enriched by flow cytometry were used as control. RNA sequencing was performed by Novogene (Nanjing, China).

We based on Cell Ranger (Version 3.0.2) pipeline coupled with a joint reference of both the mouse reference (Version mm10) and the RABV genome (NCBI record: HQ891318.1) to generate raw gene expression matrices for each sample. The Seurat 15 package (Version 3.0.0) which is embedded in R software (Version 3.5.3) was then used to analyse the output filtered gene expression matrices. First, we filtered genes expressed at a proportion <0.1% of the data and cells with <200 genes detected. Cells were identified as low-quality cells and removed from the further analyses if they met the following criteria: (1) <500 or >70 000 Unique Molecular Identifier (UMIs); (2) <200 or >7500 genes; or (3) >10% UMIs derived from mitochondrial genes. Second, we normalized the gene expression matrices and calculated 2000 features with high cell-to-cell variation by the NormalizeData and FindVariableFeatures function, respectively. Third, data was linear-transformed to scaleddata by ScaleData function, and the RunPCA function was conducted to reduce datasets dimensionality with default settings. Next, as recommended by the Seurat developers, we identified the true dimensionality of each dataset by using the ElbowPlot, DimHeatmap and JackStrawPlot functions. Finally, cells were clustered by the FindNeighbors and FindClusters functions, and the RunUMAP function was used to perform non-linear dimensional reduction with default parameters. All details regarding the scRNA-seq data processing and analysis performed in this work were adapted from the Seurat website tutorial (https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html).

The integration of multiple scRNA-seq datasets

We employed the Seurat (Version 3.0.0) integration methods to integrate multiple scRNA-seq datasets for an unbiased comparison of the changes of cell types and proportions across different conditions. Briefly, 2000 features with high cell-to-cell variation were identified as described above. Next, the FindIntegrationAnchors function was used to identify "anchors" between individual datasets and these "anchors" were inputted into the IntegrateData function to create a "batchcorrected" expression matrix of all cells. Details for the above processes refer to Seurat website tutorial (https://satijalab.org/seurat/v3.0/integration.html).

For cells identified as macrophages, all cells were first extracted from the overview integrated dataset. Next, we integrated the major cell types into an unbatched and integrated dataset for further sub-clustering. Genes were then scaled to unit variance. Finally, Scaling, PCA and clustering methods were performed for dimensional reduction as described above. The processes for sub-clustering NK cells were similar to macrophages.

Identifying cluster markers and annotating cell type After dimensional reduction, cells were gathered together in Uniform Manifold Approximation and Projection (UMAP) two-dimensional space according to their common features. We used FindAllMarkers function which is based on "wilcox" method to find markers for each of the identified clusters with parameters "min.pct = 0.2" and "logfc.threshold = 0.25". To classify and annotate different clusters, we investigated the expression of canonical markers of particular cell types. Cells expressing two or more canonical cell-type markers were identified as doublet cells and excluded from further analysis.

We performed the differential gene expression testing with the Seurat embedded function FindMarkers. To estimate the false discovery rate (FDR), the Benjamini-Hochberg method was employed. Genes with a minimum fold change of 1.5 and a maximum FDR value of 0.01 were identified as DEGs. To investigate the functions of the DEGs, we used the clusterProfiler 17 in default parameters and based on the Biological Process of GO Ontology gene sets. The Benjamini-Hochberg method was also conducted to estimate FDR in functional enrichment analysis.

Calculating cell scores to evaluate different cell states Cell scores representing the degree to which individual cells expressed a certain pre-defined expressed gene set which were initially based on the average expression of the genes. 16 To evaluate the different states of cells, AddMod-uleScore function in Seurat was used with default parameters. We used APOPTOTIC SIGNALING PATHWAY (GO: 0097190) to define the apoptosis score.

We applied the Monocle 17 (version 2) algorithm to determine the potential lineage differentiation between diverse cell populations refer to the tutorial here (http://coletrapnell-lab.github.io/monocle-release/docs/). First, store data in newCellDataSet object with the parameter "expres-sionFamily = negbinomial.size()" and "lowerDetection-Limit = 0.5". Before constructing single-cell trajectories, size factors and dispersions were estimated and filtered low-quality cells with default settings. Then the trajectory was inferred with the default parameters of Monocle after dimension reduction and cell ordering based on top 500 genes differing between clusters. Finally, the results of inferred pseudotime trajectory were presented and shown with the first two components.

The statistical tools, methods and thresholds for each analysis are explicitly described in the figure legends. No statistical methods were used to pre-determine the sample sizes. The sample size used in our study is sufficient to provide stable single-cell clustering results and to perform statistical analysis.

To map the spatial distribution of RABV in the whole brain, we constructed and characterized a recombinant RABV expressing EGFP (RABV-EGFP) in vitro and in vivo. The viral growth kinetics in BSR cells and pathogenicity in mice indicated the behaviour of RABV-EGFP was almost the same as that of the parent virus ( Figure S1 ). Then we infected C57BL/6 mice with RABV-EGFP and investigated its spread in the brain by using fMOST techniques. The research strategy was briefly depicted in Figure 1A. To evaluate the impact of infection routes on the spatial distribution of RABV in brains, groups of C57BL/6 mice (n = 10/group) were infected with different dilutions of RABV-EGFP via three different routes, including intramuscular (i.m.) injection in the hind limbs, subcutaneous (o.s.) inoculation under the ears and intranasal administration (i.n.). To be noted, i.m. injection reflects the most common infection bitten by a rabid animal in the legs; o.s. inoculation mimics the situation that the wound areas are near to the brain; i.n. administration represents several reported special cases of RABV infection caused by inhalation aerosol. 18, 19 The median lethal dose (LD 50 ) of RABV-EGFP for each infection route was calculated ( Figure S2 ). Then the mice were infected with 10×LD 50 of RABV-EGFP via three different routes. At the moribund stage, the mice were euthanized and the brains were harvested.

For fMOST analysis, the oxidized agarose-embedded brains were prepared for full-volumetric imaging with a resolution of 0.32 μm × 0.32 μm × 5 μm by automatically repeating sectioning-imaging cycles. Approximately 3000 equidistant coronal sections composed a whole-brain dataset, and the self-registration of the dataset allowed for convenient 3D visualization of the distribution of RABV-EGFP. Corresponding figures from the reference channel were used to assist the identification of brain contours and areas. To ensure the reproducibility of the technique and the veracity of the analysis, three mouse brains from each infection route were imaged. As shown in Figure 1B , we found that RABV was predominant in the cortex, hypothalamus, hippocampus, midbrain, and hindbrain. Specially, mice inoculated by i.m. infection had virus concentrated in the motor area of the cerebral cortex and the head of striatum, while mice inoculated by i.n. infection had virus primarily in dentate gyrus (DG) of hippocampal formation (HPF). The extensive RABV infection along cerebellar gyri could be observed only post i.m. inoculation. The viral load in different regions, including cerebral cortex, cerebellum, olfactory bulb, brain stem, etc., was also analysed by quantitative real-time PCR (qPCR). The results demonstrated that the viral load in olfactory bulb, brainstem, midbrain, thalamus, hippocampus and hypothalamus was varied depending on the infection route ( Figure  S3 ). Together, these findings suggest that the spatial distribution of RABV in the brain is influenced by the infection route.

The three-dimension (3D) distribution of RABV in the whole mouse brain. (A) Scheme of the experimental workflow for mapping the 3D distribution of RABV in the whole mouse brain using the fMOST technique. (B) Groups of C57BL/6 mice were inoculated with 10×LD 50 of the recombinant RABV expressing EGFP, RABV-EGFP, by the intramuscular (i.m.), the otic subcutaneous (o.s.), or the intranasal (i.n.) route. At the moribund stage, the brains were harvested and prepared for fMOST processing (n = 3). Representative pictures of the distribution of RABV in a whole mouse brain by dorsal and sagittal views are shown, and anatomical localization of the typical brain regions are shown in the left column

To further identify the brain regions infected by RABV under different inoculation routes, images of representative coronal sections, spanning from olfactory bulb to medulla, were selected and examined. The anatomical localization of RABV in these sections is shown in Table S1 . Based on cellular morphology, we found that most of RABV-infected cells were neurons, which is consistent with previous findings. 3, 20 To consolidate this finding, we stained RABV P protein (RABV-P) and neuron marker, NeuN, in the mouse brain sections, and observed them under an immunofluorescence microscope. The results showed strong co-localization between RABV-P and neurons ( Figure S4 ). To illustrate infections in a single cell, enlarged views of the motor cortex (MO), the bed nuclei of the stria terminalis (BNST) and the superior colliculus (SC) are shown in Figure 2B . The distribution of RABV in brain after i.m., o.s. and i.n. infections is shown in Figures 2C, 2D Table S1 . vestibular nuclei (VNC) and the central lobule (CENT). I.m. inoculation resulted in that RABV primarily presented in the dorsolateral part of CP, but evenly distributed throughout the CP after o.s. inoculation. In contrast, i.n. inoculation resulted in almost no RABV infection in the CP (row 1, column 5 in Figure 2C -E). RABVs were distributed along the GPe and the gyrus in the cerebellum after i.m. inoculation, but not with o.s. or i.n. routes (row 2, column 1 and 8 in Figure 2C -E).

To conduct a statistical analysis of the distribution of RABV in specific brain regions by different inoculation routes, we automatically counted the number of EGFP-positive cells in all coronal sections and calculated the cell density and number of infected cells by using the aforementioned methods (described in the Experimental Section) ( Figure 3 and Figure S5 ). To be noted, the number of EGFPpositive cells reflects the total viral load in each brain region, while the positive cell density emphasizes the concentration degree of viral infection. A majority of the cortical brain areas, except the CP and DG, showed a similar density of infection regardless of inoculation route. These areas included the paralemniscal nucleus (PL), somatosensory areas (SS), medial orbital (MO) and the infralimbic area (ILA). Similarly, the infection densities of the striatum (including accumbens nucleus (ACB), nucleus of the optic tract (OT), LS and CEA) and the hypothalamus (including PVH, AVPV, PH, lateral hypothalamic area (LHA) and ZI) were also similar among inoculation routes. In contrast, in the thalamus and brain stem, several differentially infected regions were observed, including the ventral anterior-lateral complex of the thalamus (VAL), interpeduncular nucleus (IPN), PVT, LH and VNC.

The most remarkable characteristic of rabies patients is the arousal of fear. Several nuclei in the brain have been found to be related to fear, especially CEA, BST, HPF, PAG, PVH PL, ILA, ACB, LS, SC, etc. 21 Interestingly, under all three infection routes, we identified the extensive RABV infection in these nuclei ( Figure S6B ) and the anatomical localization of these nuclei was shown in Figure S6A . Amygdala is known to play the central role in processing fear, 22, 23 and is divided into cortical division and striatal division, the latter one including medial amygdala (MEA) and CEA. Interestingly, we found that there is an extensive distribution of RABV in CEA, which is found to play a role in mediating fear-related behaviors. 24 Moreover, a newly found fear-related area, BST, 21 as well as other fear-associated regions including HPF, PAG, PVH, etc., 25 were also extensively infected by RABV. Quantification of RABV-infected cells suggested that the infection routes had no significant impact on RABV distribution in the nuclei related to fear ( Figure S6C ), indicating that these nuclei are tightly associated with the pathogenesis of rabies.

Emerging evidence has shown that RABV is not a strict neurotropic virus and besides neurons, it can infect other cell types. 26, 27 To comprehensively investigate the infection susceptibility of diverse cell populations and their contributions to RABV pathogenesis, we collected brain samples from different stages of RABV-infected mice, including mock, paralyzed and moribund mice, depending on the clinical signs post-infection. We propose that a dynamic change of viral infection in diverse cell populations will be observed. Thus, we performed droplet-based scRNAseq (10×Genomics) on a total of six mouse brain samples including two from uninfected mice (healthy: n = 2), two from mice with paralysis (paralyzed: n = 2) and two from mice in the moribund stage post RABV infection (moribund: n = 2). Following euthanasia, mouse brain tissue was obtained and rapidly digested to a single cell suspension. The single-cell suspension from the uninfected brains was directly subjected to scRNA-seq. While the single-cell suspension from RABV-infected brains was first enriched for cells containing EGFP signals (with active RABV infection) by FACS, and then analysed by scRNA-seq (Fig-ure 4A, Figure S7A ,B). With the unified single-cell analysis pipeline, ≈0.77 billion unique transcripts were obtained from 54 452 cells. Among these cells, 16 199 cells (29.75%) were from the healthy mice, 13 417 cells (24.64%) were from the paralyzed mice and 24 836 cells (45.61%) were from the moribund mice ( Figure S7C ). All high-quality cells were integrated into an unbatched and comparable dataset and subjected to principal component analysis after correction for read depth and mitochondrial read counts ( Figure  S7D ). As expected, high levels of viral mRNA in different cells were detected in both Paralyzed and Moribund conditions, and the viral load in each kind of cells is shown in Figure S7E ,F.

Using graph-based clustering, UMAP, we captured the transcriptomes of 21 high confidence cell clusters (Figure 4B) according to the expression of canonical markers ( Figure 4C) . Overall, the landscape contained the following cell lineages: T cells (Cd3d + ), natural killing (NK) cells (Klrk1 + ), B/plasma cells (Cd79a + ), monocytes (Plac8 + ), macrophages (Ms4a7 + Lyz2 + ), granulocytes (S100a8 + ), monocytes derived dendritic cells (H2-Aa + ), microglia (P2ry12 + ), neurons (Meis2 + ), oligodendrocytes (Ermn + ) and its precursors (Neu4 + ), astrocytes (Aldoc + ), ependymal cells (Tmem212 + ), Schwann cells (Ttr + ), endothelial cells (Cldn5 + ), fibroblasts (Dcn + ), meningeal cells (Slc4a10 + ) and mural cells (Vtn + ).

Consistent with previous studies, 28,29 the uninfected mouse brain consisted of microglia, neurons, endothelial, oligodendrocytes, etc. (Figure 4D ,E). We observed a great consistency in technical repeats within each sample condition, and the infected cells captured in the paralyzed and moribund conditions have shown similar patterns (Figure 4D ). After RABV infection, several cell types were significantly enriched, especially macrophages, monocytes, microglia, and NK cells ( Figure 4E ). To be noted, very few neurons were captured, possibly because neurons were very fragile post RABV infection and easy to be damaged during the processing of single-cell suspension and cell sorting by flow cytometry. To consolidate our scRNA-seq results, RABV infections in macrophages, microglia and NK cells were confirmed by confocal microscopy ( Figure  S8A ,B), virus titration ( Figure S8C ) and flow cytometry analysis ( Figure S8D ,E).

Macrophages are professional phagocytes that are integral to innate immune defence, 30 and we observed a significant enrichment of macrophages (89.3%) after RABV infection ( Figure 4E ). To gain an insight into their roles in Figure S9A ,B). Macro-C1 and Macro-C2 were predominantly found in the paralyzed and moribund mouse brains, whereas Macro-C4 was predominantly observed in the healthy mouse brains. Macro-C3 could be found in both uninfected and infected brains, while its majority resided in the paralyzed and moribund mouse brains (Figure 5B ; Figure S9A ).

To characterize each cluster, we first identified differentially expressed genes (DEGs) and then performed functional enrichment analyses ( Figure 5C,D) . Interferon response related DEGs, such as Isg15, Ifi27l2a, Ifit3, Oasl1 and Ifi209 were exclusively expressed in Macro-C1. In line with this, functional enrichment analyses revealed that Macro-C1 was an interferon response-related cluster. Macro-C2 was found to be related to cellular detoxification, metabolic and catabolic processes due to highly expressed proteases genes, such as Psmb, Psmd and Gsr. Macro-C3 highly expressed MHC-II genes, and functional enrichment analyses suggested that Macro-C3 might take part in antigen processing and presentation. Macro-C4 exclusively expressed Mrc1, Cd163, Gas6, Clec10a, Folr2 and related to myeloid leukocyte migration, similar to the recently reported border-associated macrophages (BAMs), a brain resident macrophage residing in the dura mater, subdural meninges and choroid plexus. 31 Macrophages are conventionally classified into canonical M1 and M2 classes, the pro-inflammatory and antiinflammatory macrophages, respectively. 32 We found that no macrophage cluster exhibited only M1 or M2-like phenotype, whereas Macro-C4 exhibited a more M2-dominant gene signature, such as Cd163, Mrc1 and Ccl24 ( Figure 5E ). These data indicated that macrophage activation during RABV infection did not agree with the polarization model, either as discrete states or along a spectrum of alternative polarization trajectories.

Our data probably captured infiltrating macrophages asynchronously transitioning from one transcriptomic state to the next, we thus employed Monocle2 algorithm 17 to perform the pseudotime analysis. The inferred dynamic trajectory progressing exhibited a typical branched structure: with Macro-C1 as the root, Macro-C2 and Macro-C3 as the ending clusters ( Figure 6A ; Figure S9C ). To confirm that the ordering was correct, three marker genes were selected and plotted. As shown in Figure 6B ,C, the expression level of Isg15 decreased along the pseudotime; the expression level of Psmb5 peaked in the middle of the pseudotime and the expression level of H2-DMb1 was the highest at the end of the pseudotime, demonstrating a reasonable ordering.

To understand the biological processes driving pseudotime components, we wondered which genes covary in expression with pseudotime. We clustered the representative genes identified as significantly covarying with pseudotime and identified three groups of genes expressed early, mid/mid-late and late ( Figure 6D) , consistent with the clusters identified above. Taken together, our results suggest that there are at least three different roles for the infiltrating macrophages, which we termed interferon-responsive macrophages, proteasome-active macrophages, and antigen processing and presentation macrophages during RABV infection in the brain.

NK cells are considered to be an important player of the innate immunity by controlling microbial infections. 33 Intriguingly, we found that there was an obvious enrichment of NK cells in the mouse brain after RABV infection ( Figure 4E ). Re-clustering of the total of 719 NK cells revealed that there were at least three sub-clusters (NK-C1 to C3) ( Figure 7A ; Figure S10A ,B). The NK-C1 cluster was predominantly observed in the healthy brains, whereas sub-clusters NK-C2 and NK-C3 were enriched mostly in the brains of paralyzed and moribund mice ( Figure 7B ), respectively. To be noted, the NK-C1 cluster comprised highly expressed transcription factors and genes related to gene expression regulation and chromatin assembly. The NK-C2 cluster comprised highly expressed cytotoxic genes that take part in NK cell-mediated immunity, while the NK-C3 cluster was identified as an interferon responserelated cluster ( Figure 7C,D) .

NK cell functions are regulated by a balance between activating and inhibitory signals delivered by a multitude of receptors expressed on the cell surface. 33, 34 The functional activated genes were down-regulated after RABV infection (e.g. Cd226, Klrc2, Ncr1 and Klrk1). In contrast, the exhausted genes were up-regulated after viral infection (e.g. Havcr2, Lag3, Cd96 and Ctla4). The NK cell cytotoxic genes had the highest expression levels in the NK-C2 cluster, such as Nkg7, Gzmb, Ifng and Prf1 ( Figure 7E ). The NK-C3 cluster showed the highest apoptotic degree across all three clusters ( Figure 7E,F) . These results indicate that NK cells were impaired by RABV infection with a higher exhausted and apoptotic degree than that with healthy condition.

The inferred dynamic trajectory progressing of three NK cell types also exhibited a typical branched structure with NK-C1, NK-C2 and NK-C3 distributed as the root, middle and ending clusters, respectively. We found that the pseudotime was also consistent with developmental conditions where the start corresponded to the healthy Figure S10C ). The expression level of Eomes decreased along the pseudotime, the expression level of Gzmb peaked at the middle and of Ifi205 was the highest at the end of the pseudotime (Figure 7H,I) . To understand the biological processes driving pseudotime components, we further clustered the representative genes identified as significantly covarying with pseudotime and found the marker genes clustered by pseudotemporal expression pattern also revealed the same three clusters following similar kinetic trends ( Figure 7J ).

Due to the complexity of the brain, and lacking appropriate tools, researchers are facing great difficulty to investigate the invasion and distribution of neurotropic viruses in the brain, which impedes the invention of effective therapies for neurotropic viruses. The CLARITY, approach creates a tissue-hydrogel hybrid in which tissue components are replaced with exogenous elements. And then light-microscopy techniques can be used to access the entire mouse brain. 35 Recently, 3D light sheet and confocal laser scanning microscopy was also applied to investigate the distribution of RABV in the brain and peripheral nerves. 36, 37 fMOST technology makes it possible to map the spread of neurotropic viruses in the whole brain. It rapidly acquires a full-volume dataset in a whole mouse brain with resolution to a single neuron. The recently emerging scRNA-seq provides an efficient technique to investigate the transcriptomic characteristics during infection at single-cell resolution. In this study, we investigated the spatial and cellular distribution of neurotropic virus infection in the mouse brain by jointly using fMOST and scRNA-seq. We identified specific nuclei in which viral spread or the viral load depended on the route of infection. Considering the role the infection route has on the distribution of RABV in the brain, those regions commonly infected are likely critical for deciphering RABV pathogenesis and the potential targets for future therapeutics. Considering the influence of infection route on viral distribution, knowledge of those regions that are commonly infected is especially important for understanding RABV pathogenesis. Previous studies have suggested that aggressive behaviour is critical for RABV pathogenesis, because it leads to the efficient transmission of RABV to other hosts by the rabid animals. In RABV-infected skunks, the red nucleus (RN) and midbrain raphe nuclei were found to be associated with aggressive behavior. 38 Notably, our data indicate that RABV is distributed to the RN only post i.m. inoculation ( Figure 3 ).

Fear is another cardinal characteristic of rabies, which is an emotion that has a powerful influence on behaviour and physiology. Fears, depending on the type of threat stimuli, are processed in independent neural circuits that involve the amygdala and downstream hypothalamic and brainstem circuits. 23 Rabies patients often have abnormal symptoms like hydrophobia, aerophobia and phonophobia indicating that there is a disorder in the processing of fear. In this study, we found that several nuclei related to fear, such as CEA, BNST, PAG and BLA, were extensively infected by RABV ( Figure S4) . Notably, the route of infection made no difference in the viral load in these areas, indicating that RABV infection in fear-related areas is independent of the infection route.

Besides fear, respiratory failure is another cardinal feature of rabies infection. Following rabies encephalitis, lethal complications develop quickly, including acute respiratory distress syndrome (ARDS) and myocarditis. Several cases of rabies-associated ARDS were reported, but the pathogenesis of ARDS related to RABV infection is still unknown. 39, 40 Previous studies have found that the midbrain and medulla, especially the pontine tegmentum, are the regions most affected by RABVinduced inflammation. 41 Accumulating evidence suggests that pontine nuclei, including parabrachial (PB) and Kölliker-Fuse (K-F) nuclei, contain the pontine respiratory group (PRG), which is collection of respiratory-related neurons. 42, 43 Other studies report that neurons in the SC, 45 46 also influence the respiratory cycle. In this study, we showed that the PB, SC, RN and IRN were extensively infected by RABV regardless of the infection route (Figure 3 ), suggesting that the abnormal breathing may be related to the neuronal damage caused by RABV infection in these nuclei. To be noted, RABV infection in fear and breath related regions can only provide some clue of the relationship between RABV infection and clinical signs. Further advanced techniques including optogenetics are worth to be applied to confirm this hypothesis in the future.

The fMOST results provide the spacious distribution of RABV in the brain, but could not determine the exact infected cell types. To further investigate the cellular distribution of RABV, scRNA-seq analysis was implemented and the results revealed that cell types in the brain from healthy mice consisted of microglia, neurons, endothelial, oligodendrocytes, etc. which was consistent with the previous reports. 28, 29 Different from fMOST data which showed that most of the RABV-infected cells were neurons, the vast majority of RABV-positive cells isolated from the brains of paralyzed and moribund mice were immune cells, such as macrophages and NK cells. However, the proportion of RABV-positive neurons among single-cell suspension was relatively low. A possible reason for this paradox is that neurons are hard to pass through the cell sorting by flow cytometry due to its morphology and size according to the previous studies. 46, 47 In addition, myelin and oligodendrocytes residues are very sticky and may cause neuron to cluster. Another reason is that neurons are quite fragile, especially post RABV infection, and most neurons are damaged or even die during the processing of preparing single-cell suspension, which involves repeated enzyme digestion and rotation. Before scRNA-seq analysis, we enriched EGFPpositive cells by flow cytometry and surprisingly we found that a portion of microglia /macrophages and NK cells contained RABV antigens. However, further investigation is needed to confirm that if RABV can indeed infect these cells.

Monocytes are the first line in defending against viral infection. Due to its short half-life, it is difficult for the virus to replicate in monocytes. Upon infection, monocytes change their cytokine/chemokine pattern and thus differentiate into long-lived macrophages. 48 Macrophages play a critical role in viral infection: they can restrict viral infection, while viruses can utilize macrophages as vessels for persistence or dissemination in tissues. 48 Macrophages in the brain include resident macrophages (microglia) and monocytes derived macrophages that infiltrate from the periphery post pathogen infection. CNS resident macrophages comprise microglia and BAMs. 49 We mainly observed the brain resident macrophage (Macro-C4) in healthy mice, while the macrophages in RABVinfected groups were mainly belonging to the functional macrophage types Macro-C1, Macro-C2 and Macro-C3. This is consistent with a previous study which found that peripheral monocytes quickly infiltrated into the CNS after viral infection, while resident macrophages were redundant for antigen presentation and underwent apoptosis during the chronic phase. 50 In our inferred developmental trajectory we showed that interferon responses cluster (Macro-C1) appeared earlier than cellular detoxification cluster (Macro-C2) and antigen processing and presentation cluster (Macro-C3), which is consistent with the fact that interferon induces the polarization of macrophages following virus infection, and then macrophages participate in the early immune response through phagocytosis and antigen presentation. 51 It is also noteworthy that understanding of macrophages' role in neurotropic virus infection focused on their phagocytosis and antigen presentation. However, recent studies have found that macrophages might act as a virus reservoir during infection with various viruses, assist in virus replication or long-term existence in the body, and bring viruses into other cells, even the CNS. 48 Virus reservoir macrophages are neither M1 nor M2 type, but exist in the form of M1/M2 continuum. They usually show abolished apoptosis and restricted cytopathic effects, which facili-tates the virus production. 48 In our study, we also found a group of such kind cells, suggesting that macrophages not only perform phagocytosis and antigen presentation during RABV infection, but also may play a role in assisting the virus to escape the immune response.

NK cells can eliminate infected cells by recognizing cell surface receptors or secreting killing media like perforin or TNF α . However, growing evidences showed that various viruses can induce apoptosis and lead to the exhaustion of NK cells. 52, 53 Interestingly, our inferred developmental trajectory of three NK subsets exhibited a branched structure with NK-C1, NK-C2 and NK-C3 distributed as the root, middle and ending clusters, respectively, reflecting the transition process of NK cells following virus infection from the normal state to the killing state and then to the apoptotic state. The higher apoptotic status in NK-C2 and NK-C3 than NK-C1 may indicate that RABV may escape the innate immune response by inducing the apoptosis of NK cells. In a previous study, it was found that RABV could induce the apoptosis of infiltrated T cells by upregulating FasL and B7-H1 in the surface of infected neurons, and then evading host T cell defences. 54 Thus, we speculate that RABV can employ the same strategy to evade the killing by NK cells. Recently, Huang, et al. performed scRNA-seq following intracranial injection of glycoprotein-deleted RABV to investigate the transcriptional changes in the mouse dorsal raphe nucleus (DRN). Consistent with our results, they also identified several transcriptionally diverse immune cell subsets infiltrated into the brain and uncovered their important roles after infection. 55 Further studies to investigate the variations in the types and functions of infiltrated immune cells caused by different strains of RABV with different replication capability will provide additional insights into the antiviral function of these cells.

We utilized fMOST technology to reveal the susceptible nuclei infected by RABV in the mouse brain, and take advantage of scRNA-seq technology to analyse the RABVinfected cells and illustrate the roles of some immune cells.

The joint use of these two technologies allowed us to portray an integrated map of RABV infection in the mouse brain. Our results shed a light for future investigation of the pathogenesis and clinical therapy of rabies and other neurotropic viruses such as ZIKA and Japanese encephalitis virus. 

The raw data of fMOST imaging can be accessed at http:// atlas.brainsmatics.org/a/long2111. The raw sequence data have been deposited in the Genome Sequence Archive at the BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under the accession number CRA002563 and are publicly accessible at http://bigd. big.ac.cn/gsa. Custom scripts for analyzing data are available upon reasonable request.

The authors declare no competing interests.

Bai https://orcid.org/0000-0001-6119-8936 Ling Zhao https://orcid.org/0000-0003-0569-8105

The ecology of emerging neurotropic viruses

Viral diseases of the central nervous system

The cell biology of rabies virus: using stealth to reach the brain

Micro-optical sectioning tomography to obtain a high-resolution atlas of the mouse brain

Continuously tracing brain-wide long-distance axonal projections in mice at a one-micron voxel resolution

High-throughput dual-colour precision imaging for brain-wide connectome with cytoarchitectonic landmarks at the cellular level

A platform for efficient identification of molecular phenotypes of brain-wide neural circuits

Single-cell virus sequencing of influenza infections that trigger innate immunity

Single-cell analysis identifies cellular markers of the HIV permissive cell

Single-cell profiling of ebola virus disease in vivo reveals viral and host dynamics

Potent neutralizing antibodies against SARS-CoV-2 identified by high-throughput single-cell sequencing of convalescent patients' B cells

A robust image registration interface for large volume brain atlas

A whole-brain map of long-range inputs to GABAergic interneurons in the mouse medial prefrontal cortex

A whole-brain atlas of inputs to serotonergic neurons of the dorsal and median raphe nuclei

visualization and quantification of post-stroke neural connectivity and neuroinflammation using serial two-photon tomography in the whole mouse brain

Decoding the multicellular ecosystem of lung adenocarcinoma manifested as pulmonary subsolid nodules by single-cell RNA sequencing

The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells

Effects of aerosolized rabies virus exposure on bats and mice

Airborne transmission of lyssaviruses

Rabies in the 21st century

Compensation in the neural circuitry of fear conditioning awakens learning circuits in the bed nuclei of the stria terminalis

The biology of fear

The many paths to fear

The role of the central nucleus of the amygdala in mediating fear and anxiety in the primate

Neuronal circuits for fear and anxiety

Lab-attenuated rabies virus causes abortive infection and induces cytokine expression in astrocytes by activating mitochondrial antiviral-signaling protein signaling pathway

Astrocyte infection during rabies encephalitis depends on the virus strain and infection route as demonstrated by novel quantitative 3D analysis of cell tropism

Mapping the mouse cell atlas by microwell-seq

Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris

Tissue biology perspective on macrophages

A single-cell atlas of mouse brain macrophages reveals unique transcriptional identities shaped by ontogeny and tissue environment

Much more than M1 and M2 macrophages, there are also CD169+ and TCR+ macrophages

NK cells and cancer: you can teach innate cells new tricks

Natural killer cells and other innate lymphoid cells in cancer

CLARITY for mapping the nervous system

Neuroglia infection by rabies virus after anterograde virus spread in peripheral neurons

Highresolution 3D imaging of rabies virus infection in solvent-cleared brain tissue

Diabolical effects of rabies encephalitis

Acute respiratory distress syndrome associated with rabies: a case report

Human rabies: clinical features, diagnosis, complications, and management

Pathophysiologic studies in human rabies

Functional organization of the parabrachial complex and intertrigeminal region in the control of breathing

The Kölliker-Fuse nucleus gates the postinspiratory phase of the respiratory cycle to control inspiratory off-switch and upper airway resistance in rat

Role of midbrain in the control of breathing in anuran amphibians

Hierarchy of orofacial rhythms revealed through whisking and breathing

Single-cell RNA sequencing reveals midbrain dopamine neuron diversity emerging during mouse brain development

Reduced variability of neural progenitor cells and improved purity of neuronal cultures using magnetic activated cell sorting

Monocytes and macrophages as viral targets and reservoirs

early fate defines microglia and non-parenchymal brain macrophage development

Single-cell profiling identifies myeloid cell subsets with distinct fates during neuroinflammation

Interferon-γ regulates cellular metabolism and mRNA translation to potentiate macrophage activation

Immune dysregulation in severe influenza

Severe human lower respiratory tract illness caused by respiratory syncytial virus and influenza virus is characterized by the absence of pulmonary cytotoxic lymphocyte responses

Detrimental contribution of the immuno-inhibitor B7-H1 to rabies virus encephalitis

Single-cell analysis of neuroinflammatory responses following intracranial injection of G-deleted rabies viruses