Network analysis between neuron dysfunction and neuroimmune response based on neural single-cell transcriptome of COVID-19 patients

Despite global vaccination efforts, COVID-19 breakthrough infections caused by variant virus continue to occur frequently, long-term sequelae of COVID-19 infection like neuronal dysfunction emerge as a noteworthy issue. Neuroimmune disorder induced by Inflammatory factor storm was considered as a possible reason, however, little was known about the functional factors affecting neuroimmune response to this virus. Here, using medial prefrontal cortex single cell data of COVID-19 patients, expression pattern analysis indicated that some immune-related pathway genes expressed specifically, including genes associated with T cell receptor, TNF signaling in microglia and Cytokine-cytokine receptor interaction and HIF-1 signaling pathway genes in astrocytes. Besides the well-known immune-related cell type microglia, we also observed immune-related factors like IL17D, TNFRSF1A and TLR4 expressed in Astrocytes. Based on the ligand-receptor relationship of immune-related factors, crosstalk landscape among cell clusters were analyzed. The findings indicated that astrocytes collaborated with microglia and affect excitatory neurons, participating in the process of immune response and neuronal dysfunction. Moreover, subset of astrocytes specific immune factors (hinged neuroimmune genes) were proved to correlate with Covid-19 infection and ventilator-associated pneumonia using multi-tissue RNA-seq and scRNA-seq data. Function characterization clarified that hinged neuroimmune genes were involved in activation of inflammation and hypoxia signaling pathways, which could lead to hyper-responses related neurological sequelae. Finally, a risk model was constructed and testified in RNA-seq and scRNA data of peripheral blood.


Introduction
Coronavirus disease 2019  has affected more than 432 million people (February 25th, 2022), causing more than 5.94 million people to die [1]. Besides acute damage of lung, increasing reports have shown that neurological sequelae was also frequently appeared in patients hospitalized with COVID-19(2-4). The symptoms were including loss of smell and taste, and some people may later struggle with headaches, fatigue, weakness and confusion [3,5]. It is a remarkable fact that neuron disorder could also lead to long-term neurological sequelae such as insomnia, depression and encephalitis, up to 33.62% of COVID-19 patients were diagnosed with neurological or psychiatric disorder within six months of cure [6,7]. Previously, our data also demonstrated the landscape of immune disorders with COVID-19 infected [8,9]. However, little was known about the molecular changes of neuroimmune.
It is speculated that similar with other coronaviruses like SARS and MERS, the neurological sequelae might be due to the direct penetration of COVID-19 virus on central nervous system [10][11][12]. Cytokine storms induced by infection were also considered to affect central nervous system function [13]. Severe cytokine storm has been reported in COVID-19 severity with high level of serum cytokine factors like interleukin and TNF-α(14). In infectious disease, such factors have been observed to directly across the blood-brain barrier (BBB) and active microglia and astrocytes secreting inflammatory mediators [15]. Moreover, the cytokines also compromised the BBB increasing the influx of leukocytes into brain and facilitate neuronal damage [10]. Microglia, the CNS-resident innate immune cells, Verkhratsk's review pointed out it to be in a chronic immunologically active alert state presenting low BBB protection and high vulnerability to SARS-CoV-2 infection [16]. However, it is still lack of knowledge about the landscape of microglia dysfunction involving in COVID-19 induced neuronal dysfunction, and the interactions among the cell types should be clearly charactered, especially proper molecular model for neurological sequelae prediction.
To address these questions, brain single cell data from severe COVID-19 patients were analyzed and cell clusters were identified. Through neuroimmune related ligand-receptor interaction, astrocytes were showed to cross-talk with microglia in the process of immune response, and affect excitatory neuron which contributed to neuronal disorders. The neuroimmune related hinged genes that specially expressed in astrocytes were classified and proved to correlate with Covid-19 infection and ventilator-associated pneumonia (VAP) in multi-tissues. Functional analysis clarified the hinged genes were involved in neuroimmune related pathways like activation of inflammation and hypoxia signaling pathways which further supported their significant roles in neuroimmune and neuronal dysfunction. Additionally, a risk model was conducted to predict neuronal disorders and ventilator-associated pneumonia (VAP) based on these neuroimmune hinged genes.

Data source
We obtained the latest single-cell RNA (scRNA-seq) data of medial prefrontal cortex from post-mortem data of severe COVID-19 patients (GSE159812) [17], and got the peripheral blood RNA sequencing data of COVID-19 patients from the GEO database (GSE166253), which included COVID-19 patients and healthy subjects. Next, we got the tracheal aspirate RNA-Seq sequencing data from the GEO database (GSE168017), including severe VAP COVID-19 and non-VAP COVID-19 patients. Finally, we used data of COVID-19 brain tissue (GSE164332) [18] and the scRNA-seq data from 196 COVID-19 patients with severe acute of peripheral blood and lung (GSE158055) to verify the classification efficiency of risk model. The cited information is shown in the table below. Table 1.

Single-cell data analysis of COVID-19 patients
In order to ensure high reliability and comparability of data, quality control was carried out. The process could be briefly summarized as follows [19]: 1) Filtered out cells that expressed fewer than 200 genes, because these cells might be empty droplets or formed by other abnormal causes; 2) Filtered out genes that were only expressed in less than three cells and reduced unnecessary dimensions of data; 3) Filtered the cells with significantly higher number of genes than normal value, generally filtered the top 1% of data; 4) A high proportion of mitochondria might indicate that the cell internal structure had been damaged before a single cell could build a library [20]. So filtered the cells which had more than 10% mitochondrial genes.
After the above quality control process, data was analyzed according to the standard procedure of R package "Seurat": 1) Normalized the data to eliminate the influence of library size between different cells; 2) Selected highly variable genes for analysis, generally selected the first 2000 highly variable genes. 3) PCA was used to determine the most obvious biological signals in the data set, and the first N dimension of PCA results was selected for analysis (generally n was 20); 4) Graph-based clustering method was used to divide different cells into different clusters according to cell distance matrix; 5) SFRP1(secreted frizzled related protein 1), NEUROD2(neuronal differentiation 2), GAD1(glutamate decarboxylase 1), PDGFRA(platelet derived growth factor receptor alpha), AQP4(aquaporin 4) and PTPRC(protein tyrosine phosphatase receptor type C) were used as markers to identify the major cell types in the brain: neural progenitor cells, excitatory neurons, interneurons, oligodendrocyte progenitor cells, astrocytes and microglia, respectively [19].

Functional annotation
Used the R package "clusterProfiler" to perform GO function annotation and KEGG enrichment analysis with BH-corrected p < 0.05 to select the main functional and biological pathways of genes were enriched. Similarly, we carried out functional enrichment analysis on the genes in the functional module ME-turquoise, and selected five biological pathways for display using R-package "GOplot".

Identification of differentially expressed genes
To identify genes which differentially expressed between COVID-19 patients and normal groups, we used t-test model with BH-corrected p < 0.05 and |log (FC)| > 1 to select. In addition, R-package "ggplot" was used to demonstrate the differential expression patterns of genes.

Construction of functional modules
WGCNA (WGCNA, weighted gene co-expression network analysis) could quickly extract gene co-expression modules related to sample characteristics from complex data (N multi-groups) for subsequent analysis. R packet "WGCNA" was used to construct a weight gene coexpression network. 0.85 was selected as the threshold of correlation coefficient, the optimal soft threshold was 6, and the mini mum number of module genes was selected 30 by dynamic shear tree.

Recognition of hinged genes
Based on peripheral blood data, differently expressed genes in functional module were defined as hinged genes. Cox-regression analysis was used for identifying the genes whose expression were associated with the patients' VAP. The genes with p-value for expression less than 0.05 were identified as COVID-19 VAP-related genes. Genes with hazard ratio (HR) greater than 1 were identified as risk genes and those with HR less than 1 were identified as protective genes.

Table 1
Data reference information.

Data
Data Sources Brain single cell RNA-seq GSE159812 peripheral blood RNA-seq GSE166253 tracheal aspirate RNA-Seq GSE168017 Brain RNA-seq GSE164332 blood and lung single cell RNA-seq GSE158055 2.7. Multi-factorial regulatory networks and functional characterization for hinged genes The R package "RCircos" was used to characterize the positions of hinged genes. The STRING online database (https://string-db.org/) was used to analyze the interactions among hinged genes, and the relationships with an interaction score greater than 0.7 were screened and displayed in the circos diagram. We obtained lncRNAs from lnc2target- v2.0, miRNA from miRTarBase and transcription factor (TF) targeting gene relationships from Trustv2.0 database, and used Cytoscape to construct a multi-factor regulatory network of hinged genes. Subsequently, we used GSEA software to perform GSEA enrichment analysis on the hinged genes.

Construction and evaluation of the risk model
Multivariate cox regression model was used for constructing the risk model, in which the age, sex, and race status were taken into account. Then for each patient, we calculated a risk score based on the cox regression coefficient and the expression of the gene in the multivariate cox result with p-value less than 0.05.
Where the β is the cox regression coefficient, n is the number of genes in the module, e ki is the expression level of gene k in patient i. Patients were divided into two groups based on the risk scores and log-rank test was used to evaluate the VAP difference between two groups.

Statistics and visualization of networks
All the statistical analyses were performed using R Statistical Software(V-4.0), and biological networks were visualized by Cytoscape.

Identification cell clusters in medial prefrontal cortex of COVID-19 patients
We used scRNA-seq data of medial prefrontal cortex of COVID-19 patients to investigate the molecular pathology features. 19,481 genes expressed in 5566 cells (Supplementary Figure S1) were generated by using single RNA-seq and integrated quality control pipelines (Methods), and top 5 genes of each cell were classified to get the cluster heat map (Supplementary Figure S2-A). After employing the known and universal markers [19] to authenticate neural cell populations (Supplementary Figure S2-B), six cell clusters were annotated including excitatory neurons, astrocytes, microglia, neural progenitor cells, interneurons and oligodendrocyte progenitor cells( Fig. 1-A). The proportion of excitatory neurons was the largest with 77.67%, astrocytes with 4.08% and microglia with 2.96%( Fig. 1-B). The top 20 genes in each cell cluster were mainly involved in regulation of trans-synaptic signaling, neurotransmitter transport, regulation of blood circulation, and multicellular organismal signaling. Among which we were more interested in astrocytes that were related to regulate blood circulation and multicellular organismic signaling (Supplementary Figure S3-A). In addition, some lncRNAs existed in the top 10 characteristic genes of some cell clusters, for instance EMX2OS and NEAT1 (Supplementary Figure S3-B). And it was found that lncRNAs accounted for a certain proportion in each cell cluster (Supplementary Figure S3-C), including microglia and oligodendrocyte, whose number of lncRNAs exceeded 300.
Previous studies had suggested that immune hyperactive inflammatory response was typical character of COVID-19, and cytokines storm were the main cause of multiple organ dysfunction syndrome. Interestingly, IL17D, TNFRSF1A and TLR4 expressed specifically in astrocytes ( Fig. 1-C). Therefore, to further gain the functional view of these specially expressed cytokines among different cell clusters, GO analysis showed that their related receptors were associated with viral infection response signaling pathways. There were T cell, TNF and MAPK signaling pathway in Microglia, Sphingolipid and PI3K-Akt signaling pathway in Excitatory Neurons and Neural Progenitor Cells, Cytokinecytokine receptor interaction and HIF-1 signaling pathway in Astrocytes as showed in Fig. 1-D. Notably, Hypoxia-inducible factor could regulate transcription factor of innate immunity.
According to the results, 6 main cell clusters were divided from medial prefrontal cortex of patients with COVID-19, of which excitatory neurons account for the largest proportion. Specific expression of immune-related factors was observed in different cell clusters that were functionally linked to each other, which suggests the interdependence of cell types in neuroimmune response. Interestingly, in addition to microglia as known immune cells, there were large number of immunerelated factors in astrocytes, which might be a crucial cell type for neuroimmune regulation.

Molecule characterization of Neuroimmune involving in neuronal dysfunction
Due to the frequent signal transmissions among brain cell types, we explored the interaction events between ligands and receptors from different cell types, especially focus on excitatory neurons. The results indicated that there were wide crosstalks among the main excitatory neurons, astrocytes and microglia cells (Fig. 2-A). To reveal the molecular characteristics of excitatory neurons after COVID-19 infection, we analyzed the specific expression of receptors in excitatory neurons and the corresponding ligands in astrocytes and microglia. The data showed that there was a strong crosstalk between receptor in excited neurons and ligand in astrocytes as well as microglia ( Fig. 2-B). Among them, the cytokines related genes ANGPTL4, FN1 and TGFB2 acted as ligands to affect the receptors in the excitatory neurons. We found that Hypoxia Inducible Factor 1 Subunit Alpha, HIF1A and hypoxia-related factor ANGPTL4 was highly expressed in astrocytes, and ANGPTL4 receptor immune factors ITGB1 and CDH11 were widely present in neurons ( Fig. 2-C). PRR5 as adaptive immune factor interacted with AKT1 and SGK1 is related to apoptosis, DNA repair and regulation of cell cycle. It meant that the increasing of PRR5 might promote interaction with AKT1 and SGK1, which might be factors and pathways of causing neuron disorder ( Fig. 2-D). As shown in Fig. 2-D, SGK1 was up regulated in microglia, neural progenitor and interneuron cells beside excitatory neurons, involving in multi-cell types. Fig. 2-E showed transforming growth factor-beta2 (TGFB2) highly expressed in astrocytes, which could network with ACVR1B as well as ACVR1C in other cell types. ACVR1B was widely expressed in neurons, but ACVR1C was expressed in part of excitatory neurons. And then, comparing the changes of KEGG data in excitatory neurons with ACVR1C expression or none, we found the obvious differences between the two cell sub-types (Supplementary  Table S1 and S2). It got our attention that HIF-1, Inflammatory, and Long-term depression had specifically in excitatory neurons with ACVR1C expression, as shown in Table 2. Additionally, Laminin family genes (LAMC1/LAMB2/LAMA1/LAMA4) affected the function of SV2B mainly in excitatory neurons (Fig. 2-F), which was responsible for increasing the performance of synaptic vesicles and may function in the regulation of vesicle trafficking and exocytosis [21]. While another gene regulated by LAMA1, ITGA2 was rarely distributed in excited neurons (Fig. 2-F). In short, dysregulation of the ligands specifically expressed in astrocytes affected dysfunction of the corresponding receptors in neurons, which should be crucial molecular incident for neurological sequelae. Table 2 Microglia were the first and primary immune defense cells in the central nervous system (CNS) [22]. Therefore, we analyzed the interactions between astrocytes and microglia. Our results showed that the ligands and receptors specifically expressed in astrocytes had extensive connections with microglia ( Fig. 3-A). BAG3 was specifically expressed in astrocytes, and it regulated factors HSPA9 and HSPA4 widely presented in all of cell types (Fig. 3-B). Similarly, specifically expressed ligands DLL1 (Fig. 3-C) and SPP1 (Fig. 3-D) were found in microglia cells. The factors regulated by these ligands played an important role in cell differentiation, apoptosis, and proliferation. Interestingly, these factors were rarely found in neurons, but more widely distributed in astrocytes. These results suggested that through immune related ligand-receptor interactions, astrocytes could cross-talk with microglia and functionally be involved in the process of neuroimmune disorder and neuronal dysfunction after Covid-19 infection.

Recognition of hinged Neuroimmune genes across mutli-tissues RNAseq data
According to above analysis, neuroimmune-related genes involved in neuronal dysfunction were obtained. To investigate whether these genes had responsiveness changes in other tissues of COVID-19 patients, we collected the peripheral blood RNA sequencing data and found 2324 genes were differentially expressed between COVID-19 patients and healthy subjects (Supplementary Figure S4-A). Among them, 143 genes had differential expression pattern in astrocytes cluster, involving in neuroimmune regulation and function of excitatory neurons, which could be a cluster related to neuronal -dysfunction (Supplementary Table S3). Some of these genes were confirmed to participated in the development and therapy of COVID-19, such as DPP4 was the potential treatment target to prevent COVID-19 vascular complications [23], and MARK3 affected the COVID-19 development by intervening in immune response [24]. Moreover, both lncRNA NNT-AS1(25) and PRKCQ-AS1 (26) could regulate cell proliferation, apoptosis and migration. Some differentially expressed genes such as TLR4 were considered be imperative for inflammatory response [27], which was down-regulated in peripheral blood compared to healthy individuals (Supplementary Figure S4-B) and specifically in astrocytes of COVID-19 patients ( Fig. 1-C). We also found that some genes were differentially expressed in the peripheral blood data, and they were also differentially present in different cell types of the brain (Supplementary Figure S4-C). The observations supported that the neuroimmune genes could participate in the immune response of Covid-19 infection, which further suggested their functional roles in subsequent neuronal disorders.
Meanwhile, to investigate whether such functional genes could predict COVID-19 severity, co-expression network module was constructed by the WGCNA method (Methods) based on RNA-Seq sequencing data from the tracheal aspirate of severe patients. Among them, the module ME-turquoise was prominently correlated with the severity (VAP) of COVID-19 patients (r = 0.34, p = 8E-04), so we defined ME-turquoise as the VAP developing module (Fig. 4-A). Functional enrichment analysis showed that they were significantly associated with ion transport, apoptotic signaling, and blood clotting function (Fig. 4-B). Subsequently, 14 genes were identified to correlate with VAP progressing (Methods) in Fig. 4-C. In particular, PRR5 as an immune-related factor (Fig. 4-D), sleep regulator gene-GEM (Fig. 4-D), and KCTD15 as a regulator of synaptic signal transduction (Supplementary Figure S5-F) had more significance. Next, we found that a group of genes in VAP developing module also had obvious expressional difference between patients and healthy subjects in peripheral blood transcriptome, what's more exhibited directly interaction with each other (Fig. 5-A).
Accordingly, immune factors in astrocytes could predict the development of infection and VAP in patients with peripheral blood and tracheal aspirate data. This finding supported their functional roles in COVID-19 immune response and displayed the potential involvement in neuronal injury response in the brain. Therefore, these group genes were defined as hinged neuroimmune genes of COVID-19, which should be critical factors related to neuronal dysfunction.

Multi-factorial regulatory networks and functional characterization for hinged Neuroimmune genes of COVID-19
A multi-factor regulatory network was constructed (Supplementary Figure S6-A) including transcription factors (TFs), micRNAs and lncRNAs according to method of Multi-factorial regulatory networks construction (Methods). We observed some famous regulators appearing in the hinged gene regulation network, for instance the transcription factor STAT3 acted as a transcriptional activator to affect the expression of cytokines and growth factors [28,29]. We performed GSEA (Methods) on hinged genes. Hypoxia pathway was significantly enriched (Fig. 5-B), which would cause neuronal dysfunction [30][31][32]. And the markedly enriched gene ANGPTL4 was involved in the regulation by the expression of ncRNA LAST and transcription factors including PPARA, PPARD and SMAD3 (Fig. 5-C). Besides, these genes were significantly enriched in the inflammatory related pathway-TNF-α signaling via NF-kB ( Fig. 5-D), the pathway was proved to modulate viral infection response [33]. Atypical Chemokine Receptor 3, ACKR3 in TNF-α pathway was also included in the regulation network that we found above (Fig. 5-E). In addition, these hinged genes were also enriched in IL-2/STAT5 and mTORC1 signaling pathways dramatically (Supplementary Figure S6-B, C). Overall, these results suggested that hinged neuroimmune genes and the related regulation factors were correlated with significant biological process of COVID-19.

Construction and evaluation of a neural injury risk model based on hinged neuroimmune genes
The above results suggested that the hinged neuroimmune genes were involved in important biological processes in brain, peripheral blood, and tracheal aspirate of COVID-19 patients. In particular, the severity progressing of patients could be predicted according to their expression. Hence, we constructed a risk model (Fig. 6-A) based on hinged genes (Multivariate Cox regression) [34]. The model was further confirmed to efficiently predict accumulated risk of neuronal dysregulation and VAP (Fig. 6-B). Subsequently analysis indicated that, it had significant differences in risk scores between patients younger than 60 years and those older than or equal to 60 years ( Fig. 6-C). There were also obvious differences between male and female. (Fig. 6-D). Paired survival status and ethnic again exhibited different risk scores (Fig. 6-E, F). There were eight genes in the risk model, including S100A13, PRR5, KCTD15, GPC4, GEM, ERLIN2, BAG3 and ANGPTL4 (Fig. 6-A). The genes were expressed differently in cell types of COVID-19 patients' medial prefrontal cortex, particularly in Astrocytes and Microglia cells, involving in hypoxia and inflammation response. Dysregulation of these genes affected the activity of excitatory neurons and caused cell damage response. As listed in Supplementary Table S4, the eight genes were reported to have nervous related functions and part of which were correlated with neuroimmune regulations, even directly neuron disorder. We also used RNA-seq data of COVID-19 brain tissue, getting from GSE164332 to verify the efficiency of the risk model. The risk scores were notably different between patients with COVID-19 and healthy subjects (Fig. 6-G). There were 87.5% COVID-19 patients in high-risk group and 75% normal groups in low-risk group (Fig. 6-H).
Finally, we used published single-cell data (including peripheral blood and lung tissue) to validate our model [8]. Retrieval the GSE158055 dataset, results demonstrated that the genes in our model expressed specifically in cell sub-types of innate immunity and adaptive immunity which consistent with their functions, We also found that these genes had high expression levels in COVID-19 patients with severity, but rare expression in healthy subjects ( Fig. 7 and Supplementary Fig. S7). The results indicated that through prediction model  based on hinged neuroimmune genes, it is possible to predict the neuronal dysfunction and severity of patients using expressional data from peripheral blood, lung tissue and even tracheal aspirate.

Discussion
Accumulating evidence suggests that Covid-19 infection could also cause long-term neurological sequelae among large percentage of cases along with the virus spread and variants in the world [39][40][41], besides multiorgan dysfunction and even subsequently die that caused by VAP was analyzed by using expression of PRR5 and GEM. [7,38]. The neurological sequelae is ranging from headache, anosmia and ageusia, meningitis, encephalitis to stroke [12]. This may be closely related to the immune dysregulation of specific functional brain regions, so neurological autoimmune disorders may be consequences of the SARS-CoV-2 infection. Some studies showed that cytokine storm was considered to be an important cause of COVID-19 infection and long-term sequelae and the molecular changes in immune responsiveness of the nervous system were crucial to clinical neurological dysfunction syndromes [35,36]. Therefore, detailed analysis of the neuro-immunological effects in COVID-19 infection is necessary for the prediction and treatment of neurological disorder sequelae.
So herein, we gathered and analyzed the latest GSE159812, GSE166253 and GSE168017 data within the medial prefrontal cortex, peripheral blood and tracheal aspirate. We relied on the known and universal marker genes, NEUROD2, AQP4, PTPRC, SFRP1, GAD1 and PDGFRA, to identify the major cell types in the brain: excitatory neurons, astrocytes, microglia, neural progenitor cells, interneurons and oligodendrocyte progenitor cells, respectively [19]. We found desired discrimination between different cell types. Even though the proportion of excitatory neurons was the largest, interest in astrocyte has increased dramatically because of their newly discovered roles in cytology, elongation, efficacy and plasticity [37][38][39][40]. Physiological activity of neurons triggered astrocyte signaling and the signaling from microglia and astrocytes to neurons was also sufficient to alter synapses [38]. According to ligand-receptor interaction data from scRNA-seq, we revealed that astrocytes collaborated with microglia, adjusting neuronal dysfunction through neuroimmune for the first time. The alternation of stimulated microglia could explain some of the neurological symptoms associated with COVID-19, especially fatigue, depression and 'brain fog', including confusion and forgetfulness. Interestingly, in our study results, interactions between the ligands specifically expressed in astrocytes including LAMA family genes and cytokine related genes (ANGPTL4, FN1 and TGFB2), and the network factors specifically expressed in the largest population of excitatory neurons accelerated the communication of brain neurons. The result of affections on excitatory neurons by microglia and astrocyte provided the latest evidence of neuronal dysfunction. Astrocytes also displayed extensive crosstalk with the known microglia during the neurological immune response. Our data declared that the response changes of astrocyte as well as microglia were key events for neuron complication for COVID-19.
Moreover, we successfully got a critical appraisal module of COVID-19 neurological sequelae which exhibited neuro-immunological response associated with cytokines and its receptors, using the RNAand scRNA-seq data from multi-tissue of Covid-19 patients. Comparing with scRNA data of medial prefrontal cortex, the specific expressed genes in astrocytes were associated with multicellular organismal signaling, regulation of blood circulation, trans-synaptic signaling and other functions. Marker genes in excitatory neurons were mainly related to axon genesis and regulation of cell morphogenesis. Notably, cytokines and related genes were annotated to each cluster and performed functional annotation as marker genes. These insights suggested that the neurological complications might be induced by hyper-response of cytokines network. And then we verified that cytokines and related genes had coherence of changing trend in peripheral blood. Our analysis also demonstrated that these genes were associated with hypoxic alterations, TGFB2, TNFA signaling via NF-KB and inflammation-related pathways. It had been known that neurons dysfunction usually appeared after about 4 min with hypoxia. And symptoms such as restlessness and absentminded-consciousness also would appear, which were consistent with the neurological sequelae of patients with COVID-19 [31]. Pro-inflammatory cytokine TNF had been shown to be associated with increased COVID-19 mortality [41]. Our data provided newly direct evidence that the neurological sequelae of COVID-19 were relevant to hyper-response of cytokines network.
Through single cell sequencing data analysis, we identified a group of hinged genes which encode a variety of functional factors of neuroimmunity. And the connections between changes in neuroimmune related factors and neurological disorders were established. Using such factors and their alteration pattern, a risk prediction model of neurological sequelae was built up and verified in multi-tissue data of Covid-19 patients. It also was indicated that similar changes of those factors were appeared in other infected tissues like peripheral blood and respiratory tract. So that, based on our prediction model, such hinged genes could be used in the evaluating of immune responses of nervous system, as well as in the predicting of neurological sequelae caused by Covid-19 infection.
In addition, we explored whether the genes of neuronal-dysfunction cluster were deterministic factors of severity and neurological sequelae. The constructed risk assessment model included optimal cooperatively eight genes involving in function of nervous system and signal transduction apart from neuroimmune. Such as, KCTD15 regulated neural cell formation by affecting Wnt signaling and the activity of transcription factor AP-2 [42], hypoxia-related factor ANGPTL4 acted an apoptosis survival factor [43][44][45][46], and GEM might be implicated in the proliferation of microglia and astrocytes after spinal cord injury [47]. The risk model was confirmed to efficiently predict accumulated risk of neuronal dysfunction and development of COVID-19, and indicated significant differences in risk scores by age, sex, and race. Finally, we assessed and validated this risk module by new brain tissue RNA-seq data and the latest scRNA data of peripheral blood. It was found that risk scores were different remarkably between the patients and healthy subjects, and the proportion of disease and normal samples was significantly different between the high and low risk groups.
In conclusion, our research found the crosstalking among astrocyte, microglia and excitatory neurons and revealed the molecular characterization of the neuronal dysfunction after COVID-19 infection by neural transcriptome analysis. By which, astrocyte and microglia participated in COVID-19 induced neuronal dysfunction response via hypoxia and inflammation pathways. The enrolled functional factors of astrocyte were recognized by multi-tissues transcriptome. Based on which, a molecular model was constructed that exhibited efficiently in risk prediction of neurological sequelae progressing in COVID-19 (H) Proportion differences at high-risk and low-risk group was different between COVID-19 patients (n = 9) and normal groups (n = 7), and 87.5% COVID-19 patients were in the high-risk group, in contrast, 75% normal groups were in the low-risk group. patients.

Ethics approval and consent to participate
Not applicable.

Fundings
This work was supported by Emergency Research Project on COVID-19 in Harbin Institute of Technology [XNAUEA5740100420, XNGFEQ5750000220]; the National Natural Science Foundation of Fig. 7. ScRNA data validated the model. ScRNA data from GSE158055, which included 171 patients with COVID-19, progression (n = 83)) and 25 healthy individuals. (A) ANGPTL4 expressed in CD8 cells, CD4 cells, B cells and other immunity cells, and expressed highly in COVID-19 patients with severity. (B) BAG3 expressed in immunity cells, and expressed highly in COVID-19 patients with severity. (C) GPC4 expressed in some immunity cells, and had rare expression in healthy subjects. (D) S100A13 expressed in immunity cells, and expressed highly in COVID-19 patients with severity.