key: cord-0018699-ov9vnezl
authors: Wang, Qing; Peng, Cheng; Yang, Min; Huang, Fengqi; Duan, Xuzhuo; Wang, Shaowen; Cheng, Huitao; Yang, Huirong; Zhao, Huihong; Qin, Qiwei
title: Single-cell RNA-seq landscape midbrain cell responses to red spotted grouper nervous necrosis virus infection
date: 2021-06-29
journal: PLoS Pathog
DOI: 10.1371/journal.ppat.1009665
sha: 81da3b512f2ba1de98d0f718b2d5d4157e8e0c54
doc_id: 18699
cord_uid: ov9vnezl

Viral nervous necrosis (VNN) is an acute and serious fish disease caused by nervous necrosis virus (NNV) which has been reported massive mortality in more than fifty teleost species worldwide. VNN causes damage of necrosis and vacuolation to central nervous system (CNS) cells in fish. It is difficult to identify the specific type of cell targeted by NNV, and to decipher the host immune response because of the functional diversity and highly complex anatomical and cellular composition of the CNS. In this study, we found that the red spotted grouper NNV (RGNNV) mainly attacked the midbrain of orange-spotted grouper (Epinephelus coioides). We conducted single-cell RNA-seq analysis of the midbrain of healthy and RGNNV-infected fish and identified 35 transcriptionally distinct cell subtypes, including 28 neuronal and 7 non-neuronal cell types. An evaluation of the subpopulations of immune cells revealed that macrophages were enriched in RGNNV-infected fish, and the transcriptional profiles of macrophages indicated an acute cytokine and inflammatory response. Unsupervised pseudotime analysis of immune cells showed that microglia transformed into M1-type activated macrophages to produce cytokines to reduce the damage to nerve tissue caused by the virus. We also found that RGNNV targeted neuronal cell types was GLU1 and GLU3, and we found that the key genes and pathways by which causes cell cytoplasmic vacuoles and autophagy significant enrichment, this may be the major route viruses cause cell death. These data provided a comprehensive transcriptional perspective of the grouper midbrain and the basis for further research on how viruses infect the teleost CNS.

Viral nervous necrosis (VNN) is an acute and serious fish disease caused by nervous necrosis virus (NNV) which has been reported massive mortality in more than fifty teleost species worldwide. VNN causes damage of necrosis and vacuolation to central nervous system (CNS) cells in fish. It is difficult to identify the specific type of cell targeted by NNV, and to decipher the host immune response because of the functional diversity and highly complex anatomical and cellular composition of the CNS. In this study, we found that the red spotted grouper NNV (RGNNV) mainly attacked the midbrain of orange-spotted grouper (Epinephelus coioides). We conducted single-cell RNA-seq analysis of the midbrain of healthy and RGNNV-infected fish and identified 35 transcriptionally distinct cell subtypes, including 28 neuronal and 7 non-neuronal cell types. An evaluation of the subpopulations of immune cells revealed that macrophages were enriched in RGNNV-infected fish, and the transcriptional profiles of macrophages indicated an acute cytokine and inflammatory response. Unsupervised pseudotime analysis of immune cells showed that microglia transformed into M1-type activated macrophages to produce cytokines to reduce the damage to nerve tissue caused by the virus. We also found that RGNNV targeted neuronal cell types was GLU1 and GLU3, and we found that the key genes and pathways by which causes cell cytoplasmic vacuoles and autophagy significant enrichment, this may be the major route viruses cause cell death. These data provided a comprehensive transcriptional perspective of the grouper midbrain and the basis for further research on how viruses infect the teleost CNS.

Author summary NNV is one of the most contagious and economically devastating viruses of marine fish, which lead to severe central nervous system (CNS) damage. Because of the functional a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Nervous necrosis viruses (NNVs) are non-enveloped positive-strand RNA viruses classified into the family Nodaviridae [1, 2] . NNVs are about 25-30 nm in diameter, icosahedral, nonenveloped, and have a bipartite positive-sense RNA genome. The genome of NNV is bipartite, consisting of two positive-sense RNA molecules (RNA1 and RNA2). RNA1 encodes the RNAdependent RNA polymerase (RdRp), which is responsible for viral genome replication, and RNA2 encodes the capsid protein (CP), which is the sole structural protein of NNV [3] . Fish infected with this virus show clinical signs that include abnormal swimming behavior and darkening of the fish [4] . NNVs can cause massive mortality of the larval and juvenile populations of more than 50 marine and freshwater teleost species [5] [6] [7] [8] [9] , which illustrates their strong infectivity of a wide range of hosts.

Many viruses cause serious damage to the nervous system, including Japanese encephalitis virus [10], pseudorabies virus [11], Zika virus [12], herpes simplex virus (HSV) [13] , porcine hemagglutinating encephalomyelitis virus (PHEV) [14] , and these neuronal viral infections cause high mortality. Similarly, NNV infection can lead to severe central nervous system (CNS) damage. Vacuolization of the brain can occur in infected fish and is the main cause of their death [15, 16] . Viral infection is a dynamic process driven by the interplay of antiviral cellular pathways and viral mechanisms, which have evolved to suppress antiviral activity. The membrane receptors and process of virus invasion are already well understood for many nonnervous system viruses, but little is known about the entry process and molecular signatures of CNS viruses.

The cellular heat shock cognate protein 70 (HSP70) [17] and cell surface sialic acid [18] are essential for NNV infection, and endosomal acidification also is required for effective infection [19] . However, these results do not provide information about the unique molecular signatures that arise when NNVs infect CNS cells. It is difficult to identify the specific type of cell targeted by the virus and to decipher the host immune response because of the functional diversity and highly complex anatomical and cellular composition of the CNS [20, 21] . Thus, it is important to classify the cell types present in the fish brain in order to understand the molecular mechanisms involved in viral pathogens of the CNS.

Single-cell RNA sequencing (scRNA-seq) have been applied to investigate the immune system under physiological and pathological conditions [22] [23] [24] [25] . The method provides a detailed view of the complicated immune system at single-cell resolution [26] [27] [28] . The scRNA-seq is a particularly powerful tool for identifying viral target cells, as it allows effective analysis of viral mRNAs and host signature genes in a single cell [29] [30] [31] [32] . In addition, it can provide an unbiased characterization of virus-host interactions in individual cells, which are masked at the population level [29, 30, [33] [34] [35] [36] [37] [38] .

In this study, we profiled the deep transcriptomes of tens of thousands of individual midbrain cells of orange-spotted grouper (Epinephelus coioides) harvested from control and red spotted grouper NNV (RGNNV) infected fish. We determined the transcriptional profiles of immune cells to identify those enriched in infected cells and their potential functions. Our data also indicated which cell types were targeted by RGNNV infection and the activity of specific host cell genes and pathways. The extensive datasets generated during this study will be a useful resource for detailed examination of the grouper midbrain and can serve as the basis for further research on how viruses infect the teleost CNS.

To identify areas of the brain that are damaged by RGNNV, we performed real-time quantitative polymerase chain reaction (RT-qPCR) and fluorescent in situ hybridization (FISH) to measure RGNNV content in different brain regions. The grouper (E. coioides) brain was divided into 5 parts of olfactory bulb (OB), pituitary, hypothalamus, forebrain (FB), midbrain (MB), cerebellum (CB), and spinal cord (SC). The expression levels of capsid protein (CP) and RNA-dependent RNA polymerase (RdRp) were significantly higher in the MB than in other brain regions (Fig 1A and 1B) . We generated tissue slices of the entire grouper brain and detected the replication of RGNNV in each brain region using FISH. CP-positive signals were mainly found in the MB (Fig 1C and 1D ), which indicated that this part of the brain was the main point of RGNVV infection (Fig 1D5-D7 ). The most infected areas were the torus semicircularis (TS), rostral tegmental nucleus (RT), reticular formation (RF), nucleus of the medial longitudinal fascicle (NMLF), and vascular lacuna of the area postrema (Vas) (Fig 1F-1I ).

Because we observed that RGNNV mainly infected grouper midbrain regions, we performed scRNA-seq on midbrain samples from uninfected and RGNNV-infected fish to characterize cellular heterogeneity in the midbrain (Fig 2A) . Any cell with fewer than 400 genes or more than 10% mitochondrial unique molecular identifier (UMI) counts was filtered out, and only genes with at least one UMI count detected in at least one cell were used for further analysis. After removing cells with minimum and maximum thresholds for read numbers per cell (nUMI), number of genes detected per cell (nGene), and mitochondrial RNA genes, 13,533 and 12,464 midbrain cells were obtained from one control fish (C) and one RGNNV-infected fish (NNV), respectively (S1 Table) . The nUMI, nGene, and mitochondrial mRNA percentages (pMito) are shown in S1A-S1D Fig and S1 Table. We classified cell types for all samples together and based on T-distributed stochastic neighbor embedding (tSNE) dimensionality reduction and unsupervised cell clustering. 35 cell clusters were identified based on the expressed unique transcriptional profiles as well as the top five expressed genes in each group (Fig 2B and S2 Table) .

Based on the expression of the neuronal markers synaptotagmin-1 (Syt1) and synaptosomal-associated protein 25 (Snap25), the 35 cell clusters were divided into 28 neuronal (Snap25/Syt1-high) and 7 non-neuronal clusters (Snap25/Syt1 negative or low) (Fig 2C and  2D ). The 28 neuronal clusters were further divided into 15 glutamatergic (GLU1-GLU15) and 

Landscape midbrain cell responses to RGNNV infection 13 GABAergic (GABA1-GABA18) subtypes based on the gamma-aminobutyric acid (GABA) neuronal markers Gad1 and Gad2 (Fig 2C and 2D ). Among the non-neuronal clusters, there was a subset of oligodendrocyte cells with high expression of the oligodendrocyte marker genes SRY-related HMG-box 10 (Sox10) and disrupted-in-schizophrenia 1 (Disc1) (Fig 2C  and 2D ) as well as a subset of astrocyte cells including Astro1, Astro2, and Astro3 with high expression of the astrocyte marker genes epoxide hydrolase 2 (Ephx2) and cytochrome P450 family 26 subfamily B member 1 (Cyp26b1) (Fig 2C and 2D ). There also was a subset of immune cells with high expression of the immune cell marker genes spleen focus forming virus proviral integration oncogene 1 (Spi1) and vav1 oncogene (Vav1) (Fig 2C and 2D ) and a subset of endothelial cells with high expression of the immune cell marker gene plateletderived growth factor receptor a (pdgfra) (Fig 2C and 2D ). The gene expression heat maps for the top five marker genes were generated in the 35 clusters and the genes were identified which were significantly enriched in each of the 35 clusters (S2 Fig). 

The marker genes were identified for each of the seven non-neuronal clusters. Endothelial cells were marked by glycine amidinotransferase (Gatm) and lysyl oxidase-like 3 (Loxl3) (Figs 3A and 3B and S3), and Im cells were marked by CD53 antigen (Cd53) and FYN-binding protein 1 (Fyb1) (Figs 3A and 3B and S3). The two Sox10 + cell clusters indicative of oligodendrocytes could be distinguished from one another by subtype markers myelin protein zero (Mpz) and pyridoxal phosphate binding protein (Plp) (Oligo1) and inositol polyphosphate-5-phosphatase D (Inpp5d) and apelin receptor b (Aplnrb) (Oligo2) (Figs 3A and 3B and S3). The three Cyp26b1 + astrocyte clusters could be distinguished from each other by the expression of protein patched homolog 1 (Pthc1) and PR domain zinc finger protein 16 (Prdm16) (Astro1), roundabout guidance receptor 1 (Robo1) and solute carrier organic anion transporter family member 1C1 (Slco1c1) (Astro2), and kinase binding protein (Kpb) and prosaposin receptor GPR37L1 (Gpr37l1) (Astro3) (Figs 3A and 3B and S3). The tSNE plots showed that Astro1 and Astro2 were grouped together, suggesting much more subtle differences in gene expression between these subtypes ( Fig 3B) . Astro3 formed distinct clusters, suggesting distinct molecular fingerprints ( Fig 3B) .

The expression of the non-neuronal subtype markers in grouper midbrain was confirmed by in situ hybridization (ISH) (Fig 3C) . Even though the tSNE plot showed that Astro3 (kbp) was far removed from Astro1 (pthc1) and Astro2 (Robo1), the ISH results showed that the kbp marker gene of Astro3 was expressed in the same region as Astro1 and Astro2, indicating that all three kind of cells may be the same type of cell. The Oligo 1 Plp + was mainly localized in the periventricular gray zone of the optic tectum (PDZ), RF, and NMLF areas (Fig 3C) , and the Oligo2 Aplnrb + was mainly localized in the PDZ and NMLF regions ( Fig 3C) . This indicated that Oligo1 and Oligo2 may be the same type of cell but that there were differences between them, and our marker genes distinguished between them.

Nine subpopulations of astrocytes were identified, and tSNE plots showed distinct clusters, which suggested distinct molecular fingerprints ( Fig 4A) . The gene expression heat maps were generated for the top five marker genes in nine clusters and identified the genes which were significantly enriched in each of the nine clusters ( Fig 4C and S3 Table) . Additionally, tSNE plots showed the expression of selected marker genes enriched in subsets of astrocytes ( Fig  4E) . Six subpopulations of oligodendrocytes also were identified (Fig 4B) . The gene expression heat maps were generated for the top five marker genes in 6 clusters and identified the genes which were significantly enriched in each of the six clusters ( Fig 4D and S4 Table) . The tSNE plots showed the expression of selected marker genes which were enriched in subsets of oligodendrocytes ( Fig 4F) .

The gene expression heat maps were generated and identified marker genes for each of the subpopulations of immune cells. The four subpopulations of immune cells (Figs 5A and 5B and S6 and S5 Table) could be distinguished from each other by the expression of macrophage marker gene (Mrc1) for macrophages, Pard3b for microglia, Rbfox1 for undefined, and Rhdr2 Table) , which confirmed the findings revealed by scRNA-seq. In the FISH probe test of the macrophages in the midbrain of control and RGNNV-injected fish, the Mrc1-positive signal was significantly increased after RGNNV infection ( Fig 5E) , which also was consistent with our scRNA-seq result. These findings suggested that macrophages played an important role during RGNNV infection.

To further investigate the changes in gene transcription levels of macrophages after RGNNV infection, we compared the expression patterns of genes between the control and RGNNV-infected cells. Expression of most immune-related genes was significantly up-regulated by viral infection, and these differentially expressed genes (DEGs) were involved in myeloid leukocyte activation and related to nuclear factor (NF)-κB signaling. Analysis of gene expression levels in midbrain cells from control and RGNNV-infected fish revealed that expression of most of the pro-inflammatory cytokine genes was higher in infected fish, including tumor necrosis factor alpha (TNFα), interleukin 1 (IL1), lipopolysaccharide induced TNF factor (LITAF), Toll-like receptor 18 (TLR18), CD97, and NFκB ( Fig 5E) . The RT-qPCR results also showed high expression of these genes after RGNNV infection ( Fig 5F) . Together, these findings illustrated a consistent response by innate immune cells to RGNNV infection. 

Immune cells differentiate to resist virus invasion. For example, macrophages differentiate into M1 or M2 types to exert their immune function. To evaluate this process in macrophages during RGNNV infection, we performed an unsupervised pseudotime analysis. The main cell differentiation observed was from microglia to macrophages (Fig 6A) . Furthermore, the cell trajectory branching analysis showed that Il-23, Tnf-α, and Il-1β (genes highly expressed in M1-type macrophages) were highly expressed in branch 1-2-3, but expression was significantly decreased in branch 1-2-4. Expression of Il-4, Il-10, and Tgf-β (genes highly expressed in M2 type macrophages) was significantly increased in branch 1-2-4 but significantly decreased in branch 1-2-3 ( Fig 6B and 6C) .

To further characterize the transcriptional program underlying macrophage differentiation, we then identified five groups of genes with distinct expression patterns along the differentiation process (Fig 6C) . These results suggested that RGNNV infection may lead to the transformation of microglia into macrophages, which then mainly transformed into M1-type activated macrophages.

Of the nerve cell subgroups, cell numbers in the two groups GLU1 and GLU3 were reduced the most after RGNNV infection (reduced by 40% and 47%, respectively) ( Fig 7A and 7B and S7 Table) . To confirm that RGNNV attacked GLU1 and GLU3 subgroup cells, FISH was used to test the co-localization of GLU1 and GLU3 nerve cells and RGNNV. Slc17a7, which was specifically highly expressed in both GLU1 and GLU3 subgroup cells, was used as the marker gene to mark GLU1 and GLU3 cells (Fig 7C and 7D) . The FISH results showed that Slc17a7-positive and CP-positive signals were co-located in all five RGNNV-infected fish (Fig 7E) , which indicated that RGNNV co-localized with GLU1 and GLU3 nerve cells. These results supported the premise that RGNNV mainly attacked GLU1 and GLU3 nerve cells.

To further investigate the transcriptomic changes of GLU1 and GLU3 after RGNNV infection, we compared the expression patterns of control and RGNNV-infected fish. We found that DEGs associated with membrane formation, transfer, and autophagy related pathway were significantly enriched in the RGNNV infection group, and high levels of protein tyrosine phosphatase receptor type F (ptprf), vasoactive intestinal peptide receptor 2 (vipr2), heat shock protein 90 (hsp90), lysine (K)-specific demethylase 6B (kdm6b), rho-GTPase activating protein 6 (arhgap6), and solute carrier family 6 member 6 (slc6a6) were detected (Fig 8A and 8B) . The RT-qRCR results also showed that membrane receptor, membrane lysis, transmission, and autophagy related genes were highly expressed in the RGNNV infection group (Fig 8C) .

Red spotted grouper nervous necrosis virus is a common fish pathogen that destroys the nervous system, and vacuolization of the brain is the main cause of death caused by the virus [15, 16] . Identifying the specific areas and cell types in the CNS targeted by the virus is crucial for understanding the pathogenesis of RGNNV infection. However, it is difficult to obtain an integrated scenario of the cellular and molecular immune responses that occur upon RGNNV (macrophage marker gene), and the positive signal in the RGNNV infection group was higher than that in the control group. Scale bars, 500 μm. (F) The volcano plot shows the selected DEGs from macrophages cells in the comparison of the control and RGNNV-infected fish. DEGs of control comparison with RGNNV infection group in macrophages. Each red dot denotes an individual gene with adjust P < 0.05 (P-value adjusted by false discovery rate in MAST), other indicates non-significant genes. Example genes are labeled with the gene name. (G) The gene expression levels of IFN and inflammation-related genes increased in macrophages after RGNNV infection. Data are expressed as means ± SD. ( �� p < 0.01, or ��� p < 0.001).

https://doi.org/10.1371/journal.ppat.1009665.g005

infection. To address this issue, we performed scRNA-seq and identified 35 transcriptionally distinct cell subtypes in the grouper midbrain, including 28 neuronal and seven non-neuronal cell types. Analysis of the transcriptional profiles of immune cells revealed macrophage enrichment in RGNNV-infected fish, and the macrophages were mainly polarized into M1-type macrophages and showed a strong inflammatory response. We also found that RGNNV mainly targeted GLU1 and GLU3 nerve cells and that RGNNV infection induced significant enrichment of membrane formation, lysis, and vesicular transport related pathways. These results provided a comprehensive overview of midbrain cell types based on their Researchers previously conducted a detailed census of cell types in the nervous system of mice [39, 40] , and cell types in part of the brain of reptiles also have been reported [41] . However, the types of cells in the brain of fish are poorly known, even though fish play an important role in the evolution of vertebrates. In this study, we using RGNNV infected young orangespotted grouper (E. coioides) (mean weight 3.3 ± 0.4 g, mean body length 3.7 ± 0.3 cm) and obtained 13,533 and 12,464 cells from one control fish (C) and one RGNNV-infected (NNV) fish midbrain after single cell collection. Also, we identified cell type-specific markers to allow for unambiguous cell type definition, which will enable the development of genetic or viral marker gene slc17a7 was restricted to specific GLU1 and GLU3 clusters among all of the cell clusters. (D) Violin plot showing the expression of pan marker genes across the 35 cell clusters. (E) Multicolor FISH staining to validate the Slc17a7 and RGNNV states, GLU1 and GLU3 cells were marked with Slc17a7. The middle row (1) shows high magnification of the boxed areas 1 in the top row (panorama). FISH staining shows that Slc17a7 and RGNNV were co-located in the same neuron (scale bars 10 μm). The bottom row (2) shows high magnification of the boxed areas 2 in the top row (panorama). RGNNV and Slc17A7 did not co-locate, but they localized in different areas of the same neuron (scale bars 10 μm). Scale bars in the panorama view are 200 μm.

https://doi.org/10.1371/journal.ppat.1009665.g007 tools to achieve cell type-specific labeling and manipulation. This ability is essential for analysis of cell type-specific functions in the complex midbrain region. Major classes of cells, such as neurons, astrocytes, endothelial cells, oligodendrocytes, and immune cells, can be distinguished by large sets of class-specific genes that are related to the specific function of each class of cells. In mammals and reptiles, cells in different brain regions can be divided into neuronal and non-neuronal clusters based on the expression of the pan neuronal markers Snap25 and Syt1, and neuronal clusters can be further divided into glutamatergic (GLU) and GABAergic (GABA) subtypes based on the GABA expression of Gad1 and Gad2 [41, 42] . We found the neuronal markers Snap25 and Syt1 in 28 clusters, which in turn could be divided into 15 GLU (GLU1-GLU15) and 13 GABA (GABA1-GABA13) subtypes based on their differential expression of Gad1 and Gad2. This indicated that the mammalian nervous system divisions also applied to fish [43] . We also identified potential subtype-specific marker genes for each of the 28 neuronal clusters (S2 Table) . The majority of them contained subtype-specific genes unique to that cluster, but in some cases a neuron cluster could be defined by the combinatorial expression of marker genes (Fig 2C and S2 Table) . We distinguished a number of neuron subtypes by the expression of specific genes. For example, Cbln1 and TBR1 were only found in Glu1 and Glu4 cell clusters, respectively (S2 Table) . Some genes were expressed only in certain types of nerve cells. Lhx9 was only expressed in Glu nerve cells, such as Glu6, Glu2, Glu5, and Glu13 (S2 Table) . These results demonstrated that our unbiased scRNA-seq analyses were able to identify cell types as well as cell type-specific transcriptional features in the midbrain.

In mammals, the nervous system tends to consist of both nerve cells and non-nerve cells, and non-nerve cells mainly fall into four categories according to special marker genes: astrocytes, oligodendrocytes, immune cells, and endothelial cells [40] . The types of non-nerve cells in the fish nervous system were previously unknown. Based on non-neuronal clusters marker genes in mammals and amphibians, we identified four non-nerve cell types in grouper: astrocytes (Ephx2 + and Cyp26b1 + ), oligodendrocytes (Sox10 + and Disc1 + ), immune cells (Spi1 + and Vav1 + ), and endothelial cells (Pdgfra + ) (Fig 2) [39] [40] [41] [42] . However, we also identified unique non-nerve cell marker genes in the grouper midbrain. For example, based on grouperspecific gene expression, astrocytes were divided into three clusters [Astro1 (ptch1 + and prdm16 + ), Astro2 (robo1 + and slco1c1 + ), and Astro3 (kpb + and gpr37l1 + )] and oligodendrocytes were divided into Oligo1 (plp + and mpz + ) and Oligo2 (aplnrb + and inpp5d + ). The expression of some non-neuronal subtype markers in the grouper midbrain were confirmed by ISH with riboprobes that were synthesized according the marker gene sequence (S1 Supplementary Sequence) of each non-neuronal subtype. These results suggested that the classification of CNS cells in the fish midbrain was comparable to that of other vertebrates, but the orange-spotted grouper also had its own unique marker genes.

The host immune response against acute RGNNV infection plays an antiviral role but also leads to simultaneous pathogenic injury to organs and tissues, especially in the brain of grouper. Several researchers have reported the characteristics of innate and adaptive immune responses after virus infection [43] [44] [45] [46] [47] [48] which has allowed us to speculate about the potential pathogenesis of RGNNV infection. However, it is difficult to obtain an integrated scenario of the cellular and molecular immune responses in fish infected by RGNNV. To address this issue, we profiled the immunological response landscape in RGNNV-infected grouper at single-cell resolution to identify the critical factors responsible for antiviral immunity and pathogenesis in infected fish. We identified expanded effector macrophage clusters, which may provide a possible immune response in RGNNV-infected fish. We also found that RGNNVinfected grouper showed a concerted and strong inflammatory response. Additionally, compared to control fish, viral infection led to changes in host cell transcription of several inflammatory cytokines, including TNFα, IL1, LITAF, TLR18, and NFκB. Previous studies of RGNNV showed that the virus could harness the inflammatory response to antagonize innate immunity in in vitro cell infection experiments [49] [50] [51] . In fish, interferon (IFN) is an important barrier against viral invasion [52] [53] [54] . However, none of the IFN-related genes was significantly overexpressed after RGNNV infection. Studies have shown that many RNA viruses could evade the host's immune response. For example, the rabies virus can prevent the acidification of IRF3 by TBK1 [55] , Ebola virus can inhibit the activation of IRF3 [56] , and porcine reproductive and respiratory syndrome virus suppresses type I IFN production and signaling, manipulates the cytokine responses, and modulates apoptosis to establish persistent infection in pigs [57] . A recent study in grouper spleen cells (GS) demonstrated that the B2 protein of RGNNV could inhibit the host IFN response by suppressing host transcription directed by RNA polymerase II [58] . This suggests that RGNNV may have a potential mechanism to evade the host immune response. Based on our data, we confirmed that RGNNV infection caused a concerted and strong inflammatory response in vivo, and we identified the inflammatory cytokines with high expression during RGNNV infection. However, the state of the immune system was not fully restored, as exemplified by the IFN-related subset. A long-term follow-up study is needed to investigate how fish achieve full immunity to persistent RGNNV.

In fish, macrophages are a key component of the innate immune system. Activated macrophages can produce active molecules to guide the occurrence and regression of inflammation, so they are crucial for antivirus activity and survival of organisms [59] . Microglia cells are widely considered to be immune effectors in the CNS and to play an important role in the treatment of the injured CNS. However, it is not clear whether they have any effect on the nervous system damage caused by the virus's attack [60] . Our results showed that microglia cells could be transformed into M1-type activated macrophages after viral invasion of grouper brain tissue. Studies have shown that microglia can activate macrophages to express cytokines such as IL-1 to reduce neuron damage [61, 62] . We also found that IL-1 levels were significantly increased in brain macrophages of RGNNV-infected grouper. These results suggested that viral attack of nerve cells might stimulate microglia cells to transform into M1-type activated macrophages to produce cytokines to reduce the damage to nerve tissue caused by the virus.

Previous studies reported that RGNNV mainly attacked grouper nerve cells, leading to fish death, but it was not clear which types of nerve cells are targeted by RGNNV. Previous studies also showed that the viral infection did not lead to cells being differentially excluded from the scRN-seq dataset [63] [64] [65] . Thus, in the present study we used scRNA-seq to divide nerve cells into 28 classes, and we compared control and RGNNV-infected fish to determine which cell types were lost due to virus infection. Our results showed that the number of GLU1 and GLU3 nerve cells was reduced seriously (decreases of 39.58% and 46.80%, respectively) in infected fish. Furthermore, the marker genes of GLU1/GLU3 nerve cells and virus CP genes were colocalized, which illustrated that they were the main nerve cells attacked by RGNNV. Further analysis of the genes in the GLU1/GLU3 cells that were differentially expressed between the control and RGNNV-infected fish indicated that many receptor genes were highly expressed, and these receptors may be receptors for RGNNV. Like the significantly overexpressed gene HSP90, which is a receptor for Dengue virus and Japanese encephalitis virus [66, 67] , many other receptor genes, such as Ptprf and Vipr2, were significantly highly expressed in grouper after RGNNV infection. However, identifying which receptor is the RGNNV receptor requires further study.

Previous study reported that RGNNV infection grouper spleen cells (GS) cell resulted in the fusion and enlargement of numerous cytoplasmic vacuoles [68] . However, the origin and mechanism of cytoplasmic vacuolization remain unknown. The vacuolization effects caused by viral infection have been investigated in members of 15 viral families, including hepatitis A virus, hepatitis C virus, bovine virus diarrhea virus, murine leukemia virus, Zika virus, hepatitis B virus, and polyomaviruses [69, 70] . Although the origin of virus-induced vacuoles has not been fully characterized, several reports demonstrated that the mechanisms underlying the virus-induced vacuolization effects differed among viruses. For example, 3C protease of hepatitis A virus induced numerous non-acidic cytoplasmic vacuoles, which originated from the endosome and lysosome compartments [71] . Simian virus 40 was found to induce substantial cytoplasmic vacuoles at the late productive infection stage, and binding of the viral major capsid protein VP1 to the cell surface ganglioside, GM1, triggered the formation of cytoplasmic vacuoles [72, 73] . All of these studies indicated that virus-induced vacuoles derived from different membrane organelles, including mitochondria, endoplasmic reticulum, lysosomes, Golgi apparatus, and autolysosomes. In this study, we showed that RGNNV targeted GLU1/GLU3 nerve cells, and we found that pathways related to membrane formation, lysis, and transmission were significantly enriched in these two cell types. Additionally, cytoplasmic vacuole formation and apoptosis-related genes were significantly highly expressed in these cell types after RGNNV infection. For example, lysine (K)-specific demethylase 6B (Kdm6b) was the most important gene involved in the induction of neuro-dysmorphia and it promotes cell apoptosis [74, 75] . Furthermore, rho-GTPase activating protein 6 (Arhgap6) is an important gene for cytoskeleton rearrangements and may play an important role in cytoplasmic vacuole formation [76] . Both genes were significantly highly expressed after RGNNV infection, which indicated that they may had important functions in the death of GLU1/GLU3 nerve cells. Additionally, Kaul and Lipton showed that HIV could induce injury and apoptosis in rodent and human neurons in vitro and in vivo via the release of macrophage toxic factors [77] . We found that the number of macrophages increased significantly in grouper after RGNNV infection, thus we speculated that nerve cell death may also be induced by toxic factors released by macrophages, but further research is needed.

In summary, this is the first report of the cell types present in the grouper midbrain. It provided a comprehensive transcriptional perspective of the fish midbrain, which will be a valuable resource for identifying cell type-specific functions of this complex brain region. Using scRNA-seq analysis, we demonstrated many new phenomena involved in the pathogenesis of RGNNV infection in grouper. We found that macrophages may be the key cell type required for elimination of RGNNV. Furthermore, we speculated that microglia may differentiate into M1-type activated macrophages, and this may be an immune response of grouper during RGNNV infection. Based on the analysis of DEGs of macrophages cells, we also speculated that cytokines may contribute to the persistence of RGNNV. Finally, we found that RGNNV attacked GLU1 and GLU3 cells, and we screened for genes that cause GLU1 and GLU3 cell vacuole morphology changes and death. Results of this study provided greater understanding of how fish viruses attack the nervous system and how the CNS resists viruses.

All animal experiments were carried out in strict accordance with the guidelines and regulations of the Animal Research and Ethics Committees of South China Agriculture University (SYXK-2019-0136).

Two hundred orange-spotted grouper (E. coioides) (mean weight: 3.3 ± 0.4 g, mean body length: 3.7 ± 0.3 cm) were used in this experiment. Fish were maintained at 25-30˚C in fresh seawater, which was disinfected and sterilized. The fish were fed a commercial diet according to a standard feeding scheme. For the RGNNV treatment, 50 fish were placed in each of two tanks and injected intraperitoneally with 20 μL of 107 TCID50/ml RGNNV (NNV group hereafter). For the control group (C), 50 individuals were placed in each of the two tanks and injected with 20 μL of phosphate buffered saline (PBS). when more than 50% of the fish in the virus injected group showed clinical signs of disease one control and one diseased fish were sacrificed, and midbrain were collected for single cell RNA-Seq analysis. At the same time 15 control fish and 15 diseased fish were sacrificed, and each brain was dissected and fixed in 4% buffered paraformaldehyde overnight at 4˚C, rinsed twice with cold PBS, transferred to 30% sucrose solution for 48 h, and stored at 4˚C. The samples dried with sucrose were frozen using dry ice and embedded in optimal cutting temperature (OCT) compound (Sakura, USA) and used for FISH analysis. The remaining fish (eight control and 15 showing typical disease signs) were sacrificed and whole midbrain (approximately 30 mg) tissues were frozen in liquid nitrogen immediately and stored at -80˚C before RNA extraction. This work received approval from the Animal Research and Ethics Committees of South China Agriculture University (SYXK-2019-0136).

Single-cell capture was performed using a Chromium Controller instrument (10x Genomics, Pleasanton, CA, USA), which is a highly repeatable, efficient, and stable device for cell characterization and gene expression profiling of thousands to millions of cells (https://www. 10xgenomics.com/solutions/single-cell/). Single cells were collected from the whole midbrain (approximately 30 mg) of one control fish and one RGNNV-infected fish showing typical disease signs of disease as follows. Midbrain tissue was transferred from the cryopreserved tube to a Dounce homogenizer, and 500 μL of precooled cracking buffer were added. The nuclei extraction was perform as described previously [78] . Briefly, after homogenizing the tissue, the homogenate was filtered through a 70 μm cell screen, and iodoxanol gradient solution was added prior to centrifugation. The obtained white nuclear layer was washed with nuclear cleaning buffer. After trypan blue staining, the nuclear suspension was assessed using a cell counting plate under a microscope, and the total amount, concentration, and nuclear ratio of cells with an intact nuclear membrane were calculated. The target concentration of the nuclear suspension was 700-1,200 nuclei/μL. Nucleus suspensions were loaded on the Chromium Controller instrument to generate single-cell gel bead-in-emulsions (GEMs) using Chromium Single Cell 3' Reagent v3 Kits (10× Genomics) containing a pool of~750,000 barcodes sampled to separately index the transcriptome of each cell. Thousands of individual cells were isolated into droplets together with gel beads coated with unique primers bearing 10× cell barcodes, unique molecular identifiers (UMI), and poly (dT) sequences. According to the single cell 3' reagent kit protocol, GEM-reverse transcriptions were performed in a Veriti 96-well thermal cycler (Thermo Fisher Scientific, Waltham, MA, USA). After reverse transcription, GEMs were broken and the barcoded single-strand cDNA was cleaned up using DynaBeads MyOne Silane Beads (Thermo Fisher Scientific) and a SPRI Select Reagent Kit (Beckman Coulter, Brea, CA, USA). Global amplification of cDNA was achieved using the Veriti 96-well thermal cycler, and the amplified cDNA product was cleaned up using the SPRIselect Reagent Kit.

The indexed sequencing libraries were constructed using the reagents in the Chromium Single Cell 3' Library v3 Kit for fragmentation, end repair, and A-tailing; size selection with SPRI select beads; adaptor ligation; post ligation cleanup with SPRI select beads; and sample index PCR and final cleanup with SPRI select beads. The final single cell 3' library contained the standard Illumina paired-end constructs which begin and end with P5 and P7 primers used in Illumina bridge amplification. The barcoded sequencing libraries were quantified using a Bioanalyzer Agilent 2100 System with a High Sensitivity DNA chip (Agilent, Santa Clara, CA, USA) and quantitative PCR using a KAPA Library Quantification Kit (KAPA Biosystems). Finally, two sequencing libraries were loaded onto a HiSeq2500 (Illumina, San Diego, CA, USA) with a custom paired-end sequencing mode (26 base pairs for read 1 and 98 base pairs for read 2).

The single-cell sequencing files (base call files) were processed using the Cell Ranger Single-Cell Software Suite (v2.0) for quality control, sample demultiplexing, barcode processing, and single-cell 3' gene counting [79] . The raw base call files for each sample were first multiplexed into fastq data using bcl2fastq conversion software. In the raw data analysis, the fastq data were aligned to the orange-spotted grouper genome sequence (china national genebank: CNA0000026, https://db.cngb.org/search/assembly/CNA0000026/) using STAR with default parameters. Quality control of the fastq data was performed using FastQC software, and the data were aligned to the Nucleotide Sequence Database (https://www.ncbi.nlm.nih.gov/ genbank/) using the basic local alignment search tool (BLAST) to avoid the data distortion caused by experimental contamination by other species, especially bacterial infection or contamination. After the initial quality control, the sequences with low quality barcodes and UMIs were removed.

For further counting of the UMI tags, the CellRanger count algorithm was used to generate single-cell gene counts for a single library. This approach can provide the most stable and accurate clustering solutions for 10x Genomics scRNA-seq data [80] .

Only confidently mapped, nonPCR duplicates with valid barcodes and UMIs were used to generate the gene-barcode matrix. To quantitatively identify intracellular viral segmented mRNAs to track the cells from the midbrain of fish infected with RGNNV at single-cell resolution, the scRNA-seq data of two midbrain samples were reanalyzed using the CellRanger count algorithm based the union of the mm10 and PR8 (txid211044, NCBI) reference genomes. To compare the scRNA-seq data among different libraries, the gene-cell-barcode matrix of each sample was normalized by equalizing the read depth between libraries for further merging using the CellRanger aggregate procedure, which was confirmed using the Seurat integrated analysis method [81] .

The reads from higher-depth libraries were subsampled until all libraries had an equal number of confidently mapped reads per cell. The gene-cell-barcode matrix from each of the two midbrain samples was concatenated, log-transformed, and filtered based on the number of genes detected per cell. Any cell with fewer than 400 genes or more than 10% mitochondrial UMI counts was filtered out, and only genes with at least one UMI count detected in at least one cell were used for further analysis, which was performed using CellRanger R version 2.0.0 and Seurat suite version 2.0.0.

To cluster the cells, principal component analysis (PCA) was run on the normalized filtered gene-barcode matrix to reduce the number of feature (gene) dimensions. The top five principal components were selected and passed to t-distributed Stochastic Neighbor Embedding [82] for clustering visualization in a two-dimensional space. Graph-based clustering was then run to group cells with similar expression profiles, thereby building a sparse nearest-neighbor graph without pre-specification of the number of clusters. Clusters were grouped into 35 unsupervised categories according to the differential expression profile of hallmark genes. To identify genes that were enriched in a specific cluster, the mean expression of each gene was calculated across all cells in the cluster. Each gene from the cluster then was compared to the median expression of the same gene from cells in all other clusters, and the log 2 fold-change of differentially expressed genes (DEGs) was calculated. For hierarchical clustering, pairwise Pearson's correlation between each cluster was calculated based on the mean expression of each gene across all cells in the cluster, and the log 2 fold-change of DEGs was used to create a heat map to visualize the results using MEV software (http://www.tm4.org/). Specific gene expression was graphically represented using T-distributed stochastic neighbor embedding (tSNE) plots implemented using Loupe Cell Browser software and Cell Ranger R.

The transcriptional profile data for immune cell and GLU1 and GLU3 were retrieved from the NIH SRA database with the accession code SRP040656 (https://www.ncbi.nlm.nih.gov/sra/). The most marker genes used for grouper cell classification can aligned with mouse marker gene well, most of them can identity over 60% (S8 Table) . After z-score normalization, the transcriptional profile data for immune cell and GLU1 and GLU3 from the database together with the specific cell clusters in our study were used for pairwise Pearson's correlation and PCA analysis implemented in R language (http://www.r-project.org) to demonstrate the phylogeny of specific clusters. Functional pathways representative of each gene signature was analyzed for enrichment in gene categories from the Gene Ontology Biological Processes (GO-BP) database (Gene Ontology Consortium) using DAVID Bioinformatics Resources [83] .

Monocle 2 was used to infer the pseudotemporal ordering of immune cells. We assumed the raw UMI counts were distributed according to a negative binomial distribution with fixed variance expression family to model the raw UMI count data, as recommended by the authors of Monocle 2 [83] . The Monocle 2 function BEAM was used to identify immune cells that were enriched along particular branches in the pseudotime tree. Branched heat maps were constructed using genes with q-values < 5 × 10 −5 from the BEAM.

Eight midbrains were collected from control and RGNNV-infected fish respectively. Viral RNA was extracted from samples using the QIAamp RNA Viral kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. cDNA was reverse transcribed with oligo (dT) and random hexamers using the PrimeScript RT Reagent Kit (Takara, Kyoto, Japan). RT-qPCR was performed using a Light Cycler②480 II (Roche, Basel, Switzerland). The RT-qPCR conditions were as follows: 95˚C for 30 s, followed by 40 cycles of 95˚C for 5 s, 60˚C for 30 s, and 72˚C for 30 s. Relative expression was determined by normalization to the housekeeping gene β-actin. We analyzed relative gene expression using the typical 2 -44Ct method [84] . The primers used in RT-qPCR are listed in S9 Table.

Sense and antisense digoxigenin (DIG)-labeled riboprobes were synthesized from the open reading frame sequence of the RGNNV CP, ptch1, robo1, kbp, plp, aplnrb, fyb1, and gatm genes using the DIG RNA Labeling Kit (Roche Diagnostics, Mannheim, Germany), and the biotin-labeled probes were synthesized from the open reading frame sequence of the slc17a7 gene using the biotin RNA Labeling Kit (Roche Diagnostics). The riboprobe sequences are provided in S1 Supplementary Sequence.

The procedures for RNA FISH were as follows: briefly, the brain tissues from grouper were fixed in buffered 4% paraformaldehyde overnight at 4˚C, then rinsed twice with cold PBS, transferred to 30% sucrose solution for 48 h, and stored at 4˚C. The samples dried with sucrose were frozen using dry ice and embedded in OCT compound (Sakura, USA). The tissue was sliced at -21˚C, and sections were attached to cationic anti-off slides (Thermo Fisher Scientific) and then covered and incubated at 42˚C for 1 day. Tissue sections were first prehybridized for 30 min, and then 250 μL of hybridization buffer containing 150 ng of DIG-labeled sense or antisense igf3 riboprobe were added to each slide. Slides were incubated in a humidified box at 42˚C for 16 h. After hybridization, sections were sequentially washed twice in 2×saline-sodium citrate (SSC) (1×SSC = 0.15 M NaCl, 15 mM Na citrate) at room temperature for 15 min, then in 1×SSC and 0.1×SSC at 55˚C for 1 h, and then they were mounted using Fluoroshield with DAPI (Sigma-Aldrich, St. Louis, MO, USA). Fluorescent signals from FISH were imaged using a Zeiss confocal microscope (Oberkochen, Germany).

The procedures for ISH were as follows. The slides were permeabilized three times with PBT (PBS solution containing 0.1% Tween-20) for 10 min and then treated with proteinase K in PBT (10 μg/mL) for 10 min at room temperature. They then were hybridized with sense and antisense probes at 55˚C overnight. After hybridization, sections were sequentially washed twice in 2× SSC at room temperature for 15 min and then in 1× SSC and 0.1×SSC at 55˚C for 1 h. Probes were detected with an alkaline phosphatase conjugated anti-DIG antibody (Roche Diagnostics), stained using the NBT/BCIP reagent (Roche Diagnostics), and imaged using the Zeiss microscope (Oberkochen, Germany). 

Genomic classification of fish nodaviruses by molecular phylogenetic analysis of the coat protein gene

Betanodavirus and VER Disease: A 30-year research review

Identification of host-specificity determinants in betanodaviruses by using reassortants between striped jack nervous necrosis virus and sevenband grouper nervous necrosis virus

Nodaviruses as pathogens in larval and juvenile marine finfish

Picorna-like viral particles associated with mass mortalities in larval barramundi, Lates calcarifer (Bloch)

Encephalomyelitis among turbot associated with a pocornavirus-like agent

The neurotropic herpes viruses: herpes simplex and varicella-zoster

Neurotropism of swine haemagglutinating encephalomyelitis virus (coronavirus) in mice depending upon host age and route of infection

Viral nervous necrosis in hatchery-reared larvae and juveniles of Japanese parrotfish, Oplegnathus fasciatus (Temminck & Schlegel)

A viral disease in hatchery-reared larvae and juveniles of redspotted grouper

GHSC70 is involved in the cellular entry of nervous necrosis virus

Early endocytosis pathways in SSN-1 cells infected by dragon grouper nervous necrosis virus

Inhibition of betanodavirus infection by inhibitors of endosomal acidification

A genomic atlas of mouse hypothalamic development

A new scenario of hypothalamic organization: rationale of new hypotheses introduced in the updated prosomeric model

Single-cell transcriptomics to explore the immune system in health and disease

Lung single-cell signaling interaction map reveals basophil role in macrophage imprinting

The epigenetic control of stemness in CD8+ T cell fate commitment

Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR

Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types

The neuropeptide NMU amplifies ILC2-driven allergic lung inflammation

Single-cell RNA sequencing to explore immune cell heterogeneity

Extreme heterogeneity of influenza virus infection in single cells

Dissection of influenza infection in vivo by single-cell RNA sequencing

Low ambient humidity impairs barrier function and innate resistance against influenza infection

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

The use of single-cell RNA-Seq to understand virus-host interactions

Single-cell RNA-Seq reveals transcriptional heterogeneity in latent and reactivated HIV-infected cells

Comprehensive innate immune profiling of chikungunya virus infection in pediatric cases

Defining the transcriptional landscape during cytomegalovirus latency with single-cell RNA sequencing. mBio

Single-cell transcriptional dynamics of flavivirus infection

HSV-1 single-cell analysis reveals the activation of anti-viral and developmental programs in distinct sub-populations

Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq

Molecular architecture of the mouse nervous system

Evolution of pallium, hippocampus, and cortical cell types revealed by single-cell transcriptomics in reptiles

Single-cell RNA-Seq reveals hypothalamic cell diversity

Clinical and immunological features of severe and moderate coronavirus disease 2019

Dysregulation of immune response in patients with coronavirus 2019 (COVID-19

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Innate and adaptive immunity against porcine reproductive and respiratory syndrome virus

Differential expression of innate and adaptive immune genes in the survivors of three gibel carp gynogenetic clones after herpesvirus challenge

The role of type I interferons in innate and adaptive immunity against viruses in Atlantic salmon

Red grouper nervous necrosis virus (RGNNV) induces autophagy to promote viral replication

Fish TRAF2 promotes innate immune response to RGNNV infection

Grouper PKR activation inhibits red-spotted grouper nervous necrosis virus (RGNNV) replication in infected cells

Molecular regulation of interferon antiviral response in fish

Fish interferon-stimulated genes: The antiviral effectors

Fish species-specific TRIM gene FTRCA1 negatively regulates interferon response through attenuating IRF7 transcription

Identification of the rabies virus alpha/beta interferon antagonist: phosphoprotein P interferes with phosphorylation of interferon regulatory factor 3

The Ebola virus VP35 protein inhibits activation of interferon regulatory factor 3

Regulation and evasion of antiviral immune responses by porcine reproductive and respiratory syndrome virus

ed-grouper nervous necrosis virus B2 protein negatively regulates fish interferon response by suppressing host transcription directed by RNA polymerase II

Mechanisms of fsh macrophage antimicrobial immunity

Ameboid microglia as effectors of inflammation in the central nervous system

Role of Microglial M1/M2 Polarization in relapse and remission of psychiatric disorders and diseases

Microglial activation as a compelling target for treating acute traumatic brain injury

Single-cell transcriptional landscapes reveal HIV-1-driven aberrant host gene transcription as a potential therapeutic target

Two waves of pro-inflammatory factors are released during the influenza A virus (IAV)-driven pulmonary immunopathogenesis

Early precursors and molecular determinants of tissue-resident memory CD8+ T lymphocytes revealed by single-cell RNA sequencing

Dengue virus entry into liver (HepG2) cells is independent of hsp90 and hsp70

Heat shock protein 90β in the Vero cell membrane binds Japanese encephalitis virus

Establishment of a cell line from the brain of grouper (Epinephelus akaara) for cytotoxicity testing and virus pathogenesis

Cytoplasmic vacuolization in cell death and survival

Zika virus induces massive cytoplasmic vacuolization and paraptosis-like death in infected cells

Protease 3C of hepatitis A virus induces vacuolization of lysosomal/endosomal organelles and caspase-independent cell death

Identification of a mutation in the SV40 capsid protein VP1 that influences plaque morphology, vacuolization, and receptor usage

Major capsid protein VP1 and cell surface ganglioside GM1 triggers vacuole formation

KDM6B elicits cell apoptosis by promoting nuclear translocation of FOXO1 in non-small cell lung cancer

Genetic variants in the KDM6B gene are associated with neurodevelopmental delays and dysmorphic features

Genetic screening in C. elegans identifies rho-GTPase activating protein 6 as novel HERG regulator

Chemokines and activated macrophages in HIV gp120-induced neuronal apoptosis

A single-nuclei RNA sequencing study of Mendelian and sporadic AD in the human brain

Massively parallel digital transcriptional profiling of single cells

Comparison of clustering tools in R for medium-sized 10x Genomics single-cell RNA-sequencing data

Integrating single-cell transcriptomic data across different conditions, technologies, and species

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

Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources

Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method