SARS-CoV-2 Infected Pediatric Cerebral Cortical Neurons: Transcriptomic Analysis and Potential Role of Toll-like Receptors in Pathogenesis

Different mechanisms were proposed as responsible for COVID-19 neurological symptoms but a clear one has not been established yet. In this work we aimed to study SARS-CoV-2 capacity to infect pediatric human cortical neuronal HCN-2 cells, studying the changes in the transcriptomic profile by next generation sequencing. SARS-CoV-2 was able to replicate in HCN-2 cells, that did not express ACE2, confirmed also with Western blot, and TMPRSS2. Looking for pattern recognition receptor expression, we found the deregulation of scavenger receptors, such as SR-B1, and the downregulation of genes encoding for Nod-like receptors. On the other hand, TLR1, TLR4 and TLR6 encoding for Toll-like receptors (TLRs) were upregulated. We also found the upregulation of genes encoding for ERK, JNK, NF-κB and Caspase 8 in our transcriptomic analysis. Regarding the expression of known receptors for viral RNA, only RIG-1 showed an increased expression; downstream RIG-1, the genes encoding for TRAF3, IKKε and IRF3 were downregulated. We also found the upregulation of genes encoding for chemokines and accordingly we found an increase in cytokine/chemokine levels in the medium. According to our results, it is possible to speculate that additionally to ACE2 and TMPRSS2, also other receptors may interact with SARS-CoV-2 proteins and mediate its entry or pathogenesis in pediatric cortical neurons infected with SARS-CoV-2. In particular, TLRs signaling could be crucial for the neurological involvement related to SARS-CoV-2 infection.


Introduction
COVID-19 is a severe acute respiratory disease caused by the coronavirus SARS-CoV-2. Since the first cases of infection registered in Wuhan, China in late 2019, it spread widely all over the world becoming a global pandemic. The clinical manifestations of COVID-19 can vary going from an asymptomatic infection to mild, moderate, and ultimately severe respiratory illness with multi-organ dysfunctions that may lead to the death of the patient [1].

SARS-CoV-2 Replication in Neuronal Cells
In order to verify if SARS-CoV-2 was able to replicate in neuronal HCN-2 cells we evaluated N1 and N2 copy number at 3 different time point, namely 1, 3 and 6 days post infection. We observed a time dependent increase in N1 and N2 copy numbers, indicating that the virus was able to replicate in HCN-2 neuronal cells ( Figure 1). ate the transcriptional changes induced by SARS-CoV-2 infection in pediatric human cortical neurons (HCN-2). The aim was to elucidate potential pathways involved in SARS-CoV-2 infection of neuronal cells. The discovery of SARS-CoV-2 mechanisms of infection may help to prevent nervous system early and long term complications in children and to define better pediatric treatments and prognosis.

SARS-CoV-2 Replication in Neuronal Cells
In order to verify if SARS-CoV-2 was able to replicate in neuronal HCN-2 cells we evaluated N1 and N2 copy number at 3 different time point, namely 1, 3 and 6 days post infection. We observed a time dependent increase in N1 and N2 copy numbers, indicating that the virus was able to replicate in HCN-2 neuronal cells ( Figure 1).

SARS-CoV-2 Modulated Pathways Related to Viral and Bacterial Infections
No cut-off was chosen to define the differentially expressed genes (DEGs) but only genes with q-values lower than 0.05 were inspected. The transcriptomic analysis of noninfected (HCN2-CTR) against SARS-CoV-2-infected HCN-2 cells (HCN2-SARS-CoV-2) revealed 7315 DEGs. Among them, 3527 were more expressed in HCN2-SARS-CoV-2 while 3788 were more expressed in HCN2-CTR. We enriched our data with KEGG maps and we observed 96 pathways overrepresented that are plotted in the bubble plot in Figure 2. Notably, the "Coronavirus disease-COVID-19" was among the first 10 overrepresented maps. Then, "Shigellosis", "Protein processing in endoplasmic reticulum", "Endocytosis", "AGE-RAGE signaling pathway in diabetic complications", "Focal adhesion", "Viral carcinogenesis", "Ribosome", "Salmonella infection" and "Chronic myeloid leukemia" had the highest score based on −log(q-Value). "AGE-RAGE signaling pathway in diabetic complications", "Focal adhesion", "Viral carcinogenesis", "Ribosome", "Salmonella infection" and "Chronic myeloid leukemia" had the highest score based on −log(q-Value). The pathways are vertically sorted by the ratio of the number of DEGs observed in the map over the total amount of genes that take place there. Ideally, if all the genes in the map are deregulated, the value is 1, while no DEGs in the map is 0. The size of the circle gives information about the number of DEGs included in the pathway. From the left to the right we represented a score for each pathway obtained by −log(q-Value). Interestingly, the "Coronavirus disease-COVID-19" (hsa05171) is among the first 10 overrepresented maps. Then, "Shigellosis" (hsa05131), "Protein processing in endoplasmic reticulum" (hsa04141), "Endocytosis" (hsa04144), "AGE-RAGE signaling pathway in diabetic complications" (hsa04933), "Focal adhesion" (hsa04510), "Viral carcinogenesis" (hsa05203), "Ribosome" (hsa03010), "Salmonella infection" (hsa05132) and "Chronic myeloid leukemia" (hsa05220) have the highest score. The pathways are vertically sorted by the ratio of the number of DEGs observed in the map over the total amount of genes that take place there. Ideally, if all the genes in the map are deregulated, the value is 1, while no DEGs in the map is 0. The size of the circle gives information about the number of DEGs included in the pathway. From the left to the right we represented a score for each pathway obtained by −log(q-Value). Interestingly, the "Coronavirus disease-COVID-19" (hsa05171) is among the first 10 overrepresented maps. Then, "Shigellosis" (hsa05131), "Protein processing in endoplasmic reticulum" (hsa04141), "Endocytosis" (hsa04144), "AGE-RAGE signaling pathway in diabetic complications" (hsa04933), "Focal adhesion" (hsa04510), "Viral carcinogenesis" (hsa05203), "Ribosome" (hsa03010), "Salmonella infection" (hsa05132) and "Chronic myeloid leukemia" (hsa05220) have the highest score.

SARS-CoV-2 Modulated Pattern Recognition Receptor Expression
At first, we evaluated if HCN-2 cells expressed the receptors used by SARS-CoV-2 to enter into the cells. ACE2 and TMPRSS2 that encode, respectively, for the proteins Angiotensin-converting enzyme 2 and Transmembrane protease, serine 2 were not expressed in our transcriptomic analysis. Their role in COVID-19 is widely documented because ACE2 allows the binding of SARS-CoV-2 in the lung cells [37] and TMPRSS2 acts as protein priming for the virus [38]. For this reason, the Human Protein Atlas and the Gene Expression Atlas were used to confirm that none of them are detected in the human cerebral cortex. Moreover, we performed Western blot analysis in order to verify if HCN-2 cell line expressed ACE2 protein. We also evaluated ACE2 protein in A549 cell line and in A549 cells expressing Human ACE2 (hACE2). We observed that HCN-2 cells did not express ACE2 that was instead expressed in both A549 cells and in A549-hACE2 ( Figure 3).

SARS-CoV-2 Modulated Pattern Recognition Receptor Expression
At first, we evaluated if HCN-2 cells expressed the receptors used by SARS-CoV-2 to enter into the cells. ACE2 and TMPRSS2 that encode, respectively, for the proteins Angiotensin-converting enzyme 2 and Transmembrane protease, serine 2 were not expressed in our transcriptomic analysis. Their role in COVID-19 is widely documented because ACE2 allows the binding of SARS-CoV-2 in the lung cells [37] and TMPRSS2 acts as protein priming for the virus [38]. For this reason, the Human Protein Atlas and the Gene Expression Atlas were used to confirm that none of them are detected in the human cerebral cortex. Moreover, we performed Western blot analysis in order to verify if HCN-2 cell line expressed ACE2 protein. We also evaluated ACE2 protein in A549 cell line and in A549 cells expressing Human ACE2 (hACE2). We observed that HCN-2 cells did not express ACE2 that was instead expressed in both A549 cells and in A549-hACE2 ( Figure 3). The cells use the pattern recognition receptors (PRRs) in order to recognize pathogens and start the innate immune response. Table 1 collects the PRRs deregulated in the analysis. The membrane-bound PRRs include 13 genes among which TLR1, TLR4 and TLR6 are Toll-like receptors (TLRs) while AGER, CD163L1, CD36, MEGF10, SCARB1, SCARB2, SCARF1, SSC5D, MRC1 and OLR1 are Scavenger receptors. On the other hand, 3 cytoplasmic PRRs were deregulated. DDX58 is a Rig-like receptor (RLRs) while NLRC5 and NLRP1 belong to the Nod-like receptors (NLRs) family. In addition, we showed that HCN-2 cells expressed TLR4 by western blot (Figure 3).  The cells use the pattern recognition receptors (PRRs) in order to recognize pathogens and start the innate immune response. Table 1 collects the PRRs deregulated in the analysis. The membrane-bound PRRs include 13 genes among which TLR1, TLR4 and TLR6 are Tolllike receptors (TLRs) while AGER, CD163L1, CD36, MEGF10, SCARB1, SCARB2, SCARF1, SSC5D, MRC1 and OLR1 are Scavenger receptors. On the other hand, 3 cytoplasmic PRRs were deregulated. DDX58 is a Rig-like receptor (RLRs) while NLRC5 and NLRP1 belong to the Nod-like receptors (NLRs) family. In addition, we showed that HCN-2 cells expressed TLR4 by western blot (Figure 3).

MAP2K7
122. 40  Log fold change differences for DEGs involved in immunity obtained between HCN2-SARS-CoV-2 and HCN2-CTR. For each DEG, we included the mean counts computed after DESeq2 normalization in both conditions. Furthermore, DEGs with fold change magnitude higher than 2 were inspected. Among 860 DEGs, 396 were upregulated in HCN2-SARS-CoV-2 compared to HCN2-CTR and 464 were downregulated in HCN2-SARS-CoV-2 compared to HCN2-CTR. Thus, we performed a network analysis on STRING using DEGs with fold change magnitude higher than 2 together with DEGs selected with KEGG pathways (Table 2). We used the very strict filters high confidence score and "Experiments" and "Databases" interaction sources . C5AR1, DRD4, COL9A3, MRC1, CREB5, CHRM2, UBE2V1, PYCARD, JUP,  APBB1IP, IL18R1, GADD45G, CD36, PPARG, TAS2R10, TNFRSF1B, CDON, SSTR2,  COL11A2, CDK5R1, MMP10, IRF6, PSMA2, RCSD1, BCL2, EGR1, DUSP6, PTGS1, IRF7,  PLCG2, SDC1, APO, HERC5, OLR1, NPY2R, IFI27, HTR1B, BTG2, FGF9, SALL4, REM1 Table 2. Additionally, we put a focus on the "Stress" parameter obtained after the network analysis on Cytoscape. This value gives information about the centrality of the nodes. Specifically, the stress of a node is computed counting the number of shortest paths in which the node is included. For this reason, a node with a high stress is highly involved in the biological process. In addition to the aforementioned genes, EGR1, IRF7, PPARG, SDC1 and CHD4 were the 5 DEGs with the highest stress that were not inspected in the pathways ( Figure 4). . The interactomic analysis was focused on DEGs that have a fold change magnitude higher than 2 (blue DEGs on the right). Nevertheless, we also included in the analysis DEGs that are already known to be associated with immunity response selected on the bases of KEGG analysis even if only some of them have higher fold change differences (green DEGs on the left). The size of the nodes is dependent on the stress parameter so the bigger the node size, the higher the stress of that node in the network. The stress parameter gives an index of the centrality of the node in the network. Among the blue DEGs, EGR1, IRF7, PPARG, SDC1 and CHD4 in diamond shape had the five highest stress levels. For this reason, they may have a key role in immunity response to SARS-CoV-2. . The interactomic analysis was focused on DEGs that have a fold change magnitude higher than 2 (blue DEGs on the right). Nevertheless, we also included in the analysis DEGs that are already known to be associated with immunity response selected on the bases of KEGG analysis even if only some of them have higher fold change differences (green DEGs on the left). The size of the nodes is dependent on the stress parameter so the bigger the node size, the higher the stress of that node in the network. The stress parameter gives an index of the centrality of the node in the network. Among the blue DEGs, EGR1, IRF7, PPARG, SDC1 and CHD4 in diamond shape had the five highest stress levels. For this reason, they may have a key role in immunity response to SARS-CoV-2.

SARS-CoV-2 Modulated Cytokines/Chemokines Levels in the Cell Culture Medium
The levels of cytokines/chemokines in cell culture supernatants were evaluated in order to verify the existence of a pro-inflammatory state in HCN-2 neurons infected with SARS-CoV-2. An increase in the levels of IL-1β, IL-1ra, IL-2, IL-4, IL-6, IL-17, IP-10, MCP-1, RANTES and Eotaxin was found. On the contrary, IL-10 level decreased in HCN2-SARS-CoV-2 ( Figure 5).

SARS-CoV-2 Modulated Cytokines/Chemokines Levels in the Cell Culture Medium
The levels of cytokines/chemokines in cell culture supernatants were evaluated in order to verify the existence of a pro-inflammatory state in HCN-2 neurons infected with SARS-CoV-2. An increase in the levels of IL-1β, IL-1ra, IL-2, IL-4, IL-6, IL-17, IP-10, MCP-1, RANTES and Eotaxin was found. On the contrary, IL-10 level decreased in HCN2-SARS-CoV-2 ( Figure 5).  (Table 3). IL-2, IL-4 and IL-10 levels showed only slight variations in their levels after infections and we did not find any significant gene expression modification for their transcripts or for their receptors. Additionally, a slight increase of IL-1β level was noted. In our transcriptomic analysis we found its downregulation, but its receptor encoded by IL1R1 increased (Table 3). The divergent results may depend on the increase of IL-1ra, that inhibits the activity of IL-1 binding its receptor. In line with the increase in Eotaxin, IL-6 and IL-17 levels we found in the transcriptomic analysis the upregulation of eotaxin-3, the signal transducer and receptor subunit, respectively (Table 3).   (Table 3). IL-2, IL-4 and IL-10 levels showed only slight variations in their levels after infections and we did not find any significant gene expression modification for their transcripts or for their receptors. Additionally, a slight increase of IL-1β level was noted. In our transcriptomic analysis we found its downregulation, but its receptor encoded by IL1R1 increased ( Table 3). The divergent results may depend on the increase of IL-1ra, that inhibits the activity of IL-1 binding its receptor. In line with the increase in Eotaxin, IL-6 and IL-17 levels we found in the transcriptomic analysis the upregulation of eotaxin-3, the signal transducer and receptor subunit, respectively (Table 3). Log fold change differences for DEGs encoding for cytokines/chemokines and their receptors between HCN2-SARS-CoV-2 and HCN2-CTR. For each DEG, we included the mean counts computed after DESeq2 normalization in both conditions.

SARS-CoV-2 Modulated Genes Associated to COVID-19 in UniProt
The UniProt website contains the genes that are to date associated with the infection of the SARS-CoV-2 and COVID-19 outbreaks. Among them, the 28 genes represented in Table 4 were deregulated in our analysis, so they have a potential role also in HCN-2 infection. Log fold change differences between HCN2-SARS-CoV-2 and HCN2-CTR for DEGs associated with SARS-CoV-2 infection in UniProt. For each DEG, we included the mean counts computed after DESeq2 normalization in both conditions.

Discussion
A prevalence of the nervous system's involvement in patients with COVID-19, ranging from 22.5 to 36.4% among different studies, has been described [30]. Neurological involvement of COVID-19 in neonates and children is still quite rare, but recent case reports document the potential neurologic involvement also in pediatric age [34,39]. Long period of SARS-CoV-2 infection was associated with psychiatric and cognitive disorders in adolescents and children [40]. The dormant persistence of SARS-CoV-2 in the CNS may lead to neurological complications [40].
In this work we aimed to study the capacity of SARS-CoV-2 to infect HCN-2 neuronal cortical cells, studying the changes in the transcriptomic profile. Indeed, it is not clear if SARS-CoV-2 can infect neuronal cells.

SARS-CoV-2 Replicated in Neuronal Cells and Modulated PRR Expression
Our results showed that SARS-CoV-2 was able to replicate in HCN-2 cells as suggested by the increase in copy number of N1 and N2 genes. Our results agree with a study that showed viral particles inside a human induced pluripotent stem cells-derived brain sphere model supporting the replication of the virus [41]. On the contrary, Ramani et al. found that SARS-CoV-2 targeted neurons of 3D human brain organoids but it did not appear to efficiently replicate [42]. In parallel, Tiwari et al. found that the infection of cerebral organoids was less efficient compared to lung organoids due to low ACE2 and TMPRSS2 expression [43]. On the contrary, other studies in brain organoids indicated that SARS-CoV-2 infected choroid plexus but not neurons [44,45].
Interestingly, Yi et al. found ACE2 expression in brain organoids, in the somas of mature neurons, but not in neural stem cells, but SARS-CoV-2 was observed in the axons which lack ACE2 [46]. In our experimental model, we did not find the expression for ACE2 and TMPRSS2 as confirmed also by Human Protein Atlas and the Gene Expression Atlas that reported no expression in the human cerebral cortex. Additionally, Western blot analysis showed no expression of ACE2 in HCN-2 cells. ACE2 expression in neurons is not clear and may depend on the brain area. ACE2 is expressed in choroid plexus and paraventricular nuclei of the thalamus, but not in the prefrontal cortex [27] and in mature or immature olfactory receptor neurons [47].
KEGG "Coronavirus disease-COVID-19" pathway was among the most overrepresented together with pathways related to immune system, viral and bacterial infections, as reported also in other studies [48,49].
It is possible to speculate that also other receptors may interact with SARS-CoV-2 proteins and mediate its entry or pathogenesis. For this reason, we looked at the expression of genes belonging to known PRR families, that take part in the innate immune response. The main families of PRRs are TLRs, NLRs, RLRs, and also Scavenger receptors [50]. We found the deregulation of different scavenger receptors (Table 1) and in particular one of them, SR-B1, was shown to facilitate SARS-CoV-2 entry into ACE2-expressing cells [51].
A downregulation of NLRs was found in our transcriptomic analysis (Table 1). This could indicate that SARS-CoV-2 infection may induce different cell-specific pathways as suggested by Tiwari et al. that reported induction of interferons, cytokines, and chemokines and inflammasome activation in lung organoids, while in neuronal cells SARS-CoV-2 activated TLR3/7, OAS2, complement system, and apoptotic genes [43]. Instead, TLRs and RIG-1 were upregulated. In particular, TLRs seem to be involved in COVID-19 pathogenesis and emerged as potential entry factors [29,52].
We also found the upregulation of DDX58, encoding for RIG-1, a known receptor for viral RNA, that led to transcriptional induction of genes encoding type I interferons [53]. However, downstream RIG-1, the genes encoding for TRAF3, IKKε (IKBKE), and IRF3 were downregulated. It was already reported that SARS-CoV-2 is able to antagonize IFN-I production and signaling [54][55][56]. We also found ISG15 upregulation, that is targeted by SARS-CoV-2 protease, to attenuate type I interferon responses [57]. Moreover, a negative regulation of RIG-I-mediated antiviral activity due to ISG15 conjugation was reported [58]. In our transcriptomic analysis interferon response seemed to be attenuated, in line with previous results in which human brain organoids were infected but no type I interferon response was detected [26].

SARS-CoV-2 Induced TLRs, NF-κB and MAPK Signaling
In our transcriptomic analysis we found the upregulation of TLR1, TLR4 and TLR6. The expression of TLR4 was also confirmed by Western blot analysis. Interestingly, a molecular docking study indicated the binding of native spike protein of SARS-CoV-2 to TLR1, TLR4, and TLR6 and specifically, TLR4 possessed the strongest binding affinity to spike protein [59]. Other evidence suggested that SARS-CoV-2 S-protein can bind with TLR4/MD-2 complex [60] and can activate TLR4 in monocyte, neutrophil and macrophage cell lines [61]. It is important to notice that activation of TLR4 by spike protein was not regulated by ACE2 and TMPRSS2 or virus entry [61]. These results may indicate a role for TLRs in the manifestations of COVID-19 pathology as suggested by an enhancement of TLR4 mediated inflammatory signaling in PBMCs of COVID-19 patients [62].
The canonical pathway of TLR signaling leads to NF-κB activation, that in turn increases the expression of cytokines and chemokines. In particular, TLRs, through an activation cascade, activate TAK1 and TAB2, encoded by MAP3K7 and TAB2, respectively. Activated TAK1 leads to the activation of IκB kinase (IKK), that phosphorylates and degrades the inhibitory molecule IκB [63] allowing dimerization and nuclear translocation of the p65 and p50 subunits of NF-κB. TAK1 can also activate the ERK1/2 mitogen-activated protein kinases (MAPK), JNK, and p38 [63]. Our transcriptomic analysis seemed to indicate the activation of TLR4 signaling and of its downstream mediators NF-κB and MAPK. In our transcriptomic analysis we found the upregulation of the genes encoding for TAB2 and TAK1. Moreover, we found the upregulation of CHUK encoding for the catalytic subunit IKKα, while IKBKB and IKBKG encoding for IKKβ and IKKγ were downregulated ( Table 2). We also found the upregulation of NF-κB p65 subunit (RELA) and of its inhibitor IκBα (NFKBIA). The upregulation of NFKBIA is not surprising because it depends on a negative feedback loop.
In addition, ERK (MAPK3) and JNK (MAPK10) were upregulated in our transcriptomic analysis, while p38 (MAPK12 and MAPK13) and MKK3/7 (MAP2K3 and MAP2K7) were downregulated ( Table 2). ERK activation may prevent host cell apoptosis causing viral spread and the persistence of viral infection [64]. MAP2K7 downregulation may depend on feedback mechanisms or on mutual regulation by the different pathways. Indeed, it is known that the NF-κB pathway is able to modulate JNK, acting on MKK7 [65]. The involvement of JNK instead of p38 may be a feature of infected neurons. A previous work reported that JNK cascade was activated in influenza A virus-infected neurons, while a delayed p38 activation was visible in astrocytes [66]. JUN and FOS were downregulated, as found also in a transcriptomic analysis of nasopharyngeal swabs of SARS-CoV-2 positive compared to negative subjects [67]. Interestingly, JUN transcript was restricted at post transcriptional level by SARS-CoV-2 [68], also due to its short half-life [69]. However, looking at the transcription factors activated by MAPK in KEGG, we found the upregulation of ATF4, reported to be activated by TLR4 pro-inflammatory signaling [70]. Interestingly, ATF4 and NF-κB were suggested as drivers of the proinflammatory cytokine response in human airway epithelial cells infected by SARS-CoV-2 [71].
These results indicated that SARS-CoV-2 infection of cortical neurons led to TLR activation, with the induction of NF-κB and MAPK signaling. The ability of SARS-CoV-2 S1 protein to induce pro-inflammatory mediators through NF-κB and JNK activation via TLR4 was also demonstrated in macrophages [72]. These results should also be confirmed in other neuronal cell lines, in order to verify if this mechanism is valid for other cortical neurons and for other neuronal cell types.
In the cell culture medium, we also found a slight increase in the levels of IL-1β. However, we found IL1B downregulation, but the expression of its receptor IL1R1 was upregulated. The low levels of IL-1β may depend on the increase of IL-1ra, that binds to IL1R1 and inhibits the activity of IL-1. Additionally, an increase in IL-6 and IL-17 levels in the medium was found and in parallel their receptor subunits showed an increased expression.
We found the upregulation of the gene encoding for caspase 8, that mediates cell death and inflammation caused by SARS-CoV-2 infection of lung epithelial cells [79].

Network Analysis and UniProt Repository Inspection
Then, we performed a network analysis to find the genes with fold change >2 of our transcriptomic analysis that interact directly with the selected genes. In particular, we evidenced that among the interacting genes, EGR1, IRF7, PPARG, SDC1 and CHD4 were the five with a higher stress value, that indicate their relevance in connecting the regulatory molecules. These genes were all downregulated and in particular EGR1 and IRF7 are involved in interferon signaling, PPARG and SDC1 may be involved in SARS-CoV-2 pathogenesis [80,81], while CHD4 belonged to the pathway viral carcinogenesis reported to be significant in our transcriptomic analysis.
Moreover, in order to deepen the knowledge about the SAR-CoV-2 infection of neurons, we searched for genes known to be involved in COVID-19 reviewed by UniProt. In total, 28 genes showed a significantly changed expression and some of them were involved in interferon, cytokine, and RIG-1 signaling. We also found the upregulation of IFNAR1, a receptor for Type I interferon with low affinity. This may indicate also the lack of interferon, indeed, it was reported that interferon stimulation led to the downregulation of its receptor [82]. Other genes were reported to play a role in virus entry, such as FURIN, needed for proteolytic activation of SARS-CoV-2 [83], BSG encoding for CD147, with contrasting results about its role as an entry factor [84,85], and PIKFYVE [86]. This result indicated that also in neurons, even in the absence of ACE2, SARS-CoV-2 is able to induce transcriptional changes of genes associated to its pathogenesis.
However, it is important to notice that no cut-off was chosen to define DEGs but only genes which q-value was lower than 0.05 were inspected. For this reason, some of the genes present a low, even if significant, fold change. Indeed, even if a gene shows a low fold change, it can play an important role in a biological signaling.

In Vitro HCN-2 SARS-CoV-2 Infection Assay
HCN-2 cells (ATCC CRL-10742) were cultured in DMEM (Euroclone, Milan, Italy) with 10% FBS medium and 100 U/mL penicillin and 100 µg/mL streptomycin, in a 25 cm 2 culture flask. DMEM containing 100 U/mL penicillin and 100 µg/mL streptomycin was used as inoculum in the mock-infected cells. Cell cultures were infected with 1 multiplicity of infection (MOI) and incubated at 37 • C and 5% CO 2 . After 3 h cells were washed two times with lukewarm PBS and refilled with the proper growth medium (10% FBS). Optical microscope observation (ZOE™ Fluorescent Cell Imager, Bio-Rad, Hercules, CA, USA) was performed daily to investigate the cytopathic effect. SARS-CoV-2 RNA was extracted from 1, 3 and 6 days post infection (dpi) supernatant using the Maxwell ® RSC Instrument with Maxwell ® RSC Viral Total Nucleic Acid Purification Kit (Promega, Fitchburg, WI, USA). Viral RNA was reverse transcribed in a single-step RT-qPCR (GoTaq 1-Step RT-qPCR; Promega, Fitchburg, WI, USA) on a CFX96 instrument (Bio-Rad, Hercules, CA, USA) using primers specifically designed to target two regions of the nucleocapsid (N1 and N2) gene of SARS-CoV-2 (2019-nCoV CDC qPCR Probe Assay emergency kit; IDT, Coralville, IA, USA), together with primers for the human RNase P gene. Viral copy quantification was assessed by creating a standard curve from the quantified 2019-nCoV_N positive Plasmid Control (IDT, Coralville, IA, USA). The infected cells were harvested for mRNA collection 6 dpi. RNA was extracted from mock and infected HCN-2 with the acid guanidium thiocyanatephenol-chloroform method by using the RNAzolTM B (Tel-Test, Inc., Friendswood, TX, USA) and dissolved in RNAse-free water.

RNA-seq Analysis
The library preparation was performed following the TruSeq RNA Exome protocol (Illumina, San Diego, CA, USA) following the instruction, as previously described by Silvestro et al. [88]. The Illumina MiSeq Instrument was used to obtain the raw data from the libraries. The tool fastQC (version 0.11.5, Babraham Institute, Cambridge, UK) was used to check the quality parameters of the raw data in Fastq format. We counted 17,605,130 reads for the HCN2-CTR and 10,128,701 reads for HCN2-SARS-CoV-2. The low quality bases and the adapters were filtered out by Trimmomatic (version 0.38, Usadel Lab, Aachen, Germany) [89]. The remaining reads were aligned against the GRCh38 version of the human reference genome with Spliced Transcripts Alignment to a Reference (STAR) RNA-seq aligner (version 2.7.3a, New York, NY, USA) [90]. We checked for coverage of bam files using the tool samtools [91] with the command depth. We observed a mean coverage of 4.02624 for HCN2-CTR and of 2.71703 for HCN2-SARS-CoV-2. Then, the python package htseq-count (version 0.6.1p1, European Molecular Biology Laboratory (EMBL), Heidelberg, Germany) [92] associated the counts to each transcript and the package DESeq2 of Bioconductor [93] were used to compute the analysis of differentially expressed genes in R (version 3.6.3, R Core Team). No fold change threshold was used to discard the DEGs but the post hoc Benjamini-Hochberg procedure was adopted to adjust the p-value and discard DEGs with a q-value higher than 0.05.

In Silico Data Analysis
The Human Protein Atlas (http://www.proteinatlas.org (accessed on 11 June 2021)) [94] and the Gene Expression Atlas (https://www.ebi.ac.uk/gxa/home (accessed on 11 June 2021)) [95] databases were used to confirm the absence of ACE2 and TMPRSS2 in the human cerebral cortex. In R, we searched for the pathways overrepresented in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [96] using the enirchKEGG function of the library clusterProfiler [97]. The name of each DEG was converted in its related entrez using the ensembl package biomaRt [98]. The Benjamini-Hochberg method was adopted to reduce the false positives. Only the maps with a q-value lower than 0.05 were represented in the Bubble Plot realized with the library ggplot2 [99]. To have a clear representation of the maps we replaced the q-value with a score −log(q-value). Additionally, KEGG was used to inspect the curated biological pathways related to Coronavirus Disease-COVID-19 (hsa05171), Toll-like receptor signaling pathway (hsa04620), NF-kappa B signaling pathway (hsa04064). Then, a final interactomic analysis was carried out with STRING database [100]. The high confidence (0.700) level and "Experiments" and "Databases" interaction sources were required. Finally, CYTOSCAPE (version 3.8.2, Institute for Systems Biology, Seattle, WA, USA) [101] was used to customize the network and perform network analysis. The HUGO Gene Nomenclature Committee (HGNC) website [102] was used to characterize the PRRs in the Toll like receptors (group 948), NLR family (group 666) and Scavenger receptors (1253). Additionally, the category Rig-like receptors from the PRRDB 2.0: Pattern Recognition Receptor Database [103] was included. The COVID-19 section of the UniProt (https://covid-19.uniprot.org/ (accessed on 11 June 2021)) [104] website contains the reviewed proteins associated with SARS-CoV-2.

Cytokine and Chemokine Measurement by Multiplex Assay
Concentration of cytokines/chemokines was assessed in cell culture supernatants 6 days post infection by using immunoassays formatted on magnetic beads (Bio-Plex Pro Human Cytokine 27-plex Assay #M500KCAF0Y) (Bio-Rad, Hercules, CA, USA), according to manufacturer's protocol via Luminex 100 technology (Luminex, Austin, TX, USA). Briefly, the capture antibody-coupled beads are first incubated with antigen standards or samples for a specific time. The plate is then washed to remove unbound materials, followed by incubation with biotinylated detection antibodies. After washing away the unbound biotinylated antibodies, the beads are incubated with a reporter streptavidinphycoerythrin conjugate (SA-PE). Following removal of excess SA-PE, the beads are passed through the array reader, which measures the fluorescence of the bound SA-PE. For the targets over-range an arbitrary value of 4000 pg/mL is assigned, while 0 pg/mL is attributed to values below limit of detection.

Statistical Analysis
GraphPad Prism 6.0 software (GraphPad Software, La Jolla, CA, USA) was used to perform the statistical analysis. The multiple comparison was performed using one-way ANOVA test and the Bonferroni post hoc test. The comparison between two groups was performed using Student's t-test. A p value less than or equal to 0.05 was considered statistically significant.

Conclusions
This study indicated that SARS-CoV-2 induced a pro-inflammatory response in cortical neurons even if ACE2 receptor is not expressed. The transcriptomic analysis evidenced that different PPRs are deregulated and in particular TLR1, TLR4 and TLR6 signaling is upregulated leading to the activation of NF-κB and MAPK signaling. As a consequence, COX-2, MMP3 and chemokines, that may represent a damage signal in order to enroll immune system cells, were upregulated. Moreover, we observed a deregulation of type I interferon pathway, supporting the idea that SARS-CoV-2 attenuates this pathway to evade the immune system. These results evidenced that SARS-CoV-2 can activate a proinflammatory state also in cells not expressing ACE2.