Deregulation of ceRNA Networks in Frontal Cortex and Choroid Plexus of Brain during SARS-CoV-2 Infection Aggravates Neurological Manifestations: An Insight from Bulk and Single-Cell Transcriptomic Analyses.

Although transcriptomic studies of SARS-CoV-2-infected brains have depicted variability in gene expression, the landscape of deregulated cell-specific regulatory circuits has not been elucidated yet. Hence, bulk and single-cell RNA-seq data are analyzed to gain detailed insights. Initially, two ceRNA networks with 19 and 3 differentially expressed (DE) hub lncRNAs are reconstructed in SARS-CoV-2 infected Frontal Cortex (FC) and Choroid Plexus (CP), respectively. Functional and pathway enrichment analyses of downstream mRNAs of deregulated ceRNA axes demonstrate impairment of neurological processes. Mapping of hub lncRNA-mRNA pairs from bulk RNA-seq with snRNA-seq data has indicated that NORAD, NEAT1, and STXBP5-AS1 are downregulated across 4, 4, and 2 FC cell types, respectively. At the same time, MIRLET7BHG and MALAT1 are upregulated in excitatory neurons of FC and neurons of CP, respectively. Here, it is hypothesized that downregulation of NORAD, NEAT1, and STXBP5-AS1, and upregulation of MIRLET7BHG and MALAT1 might deregulate respectively 51, 6, and 37, and 31 and 19 mRNAs in cell types of FC and CP. Afterward, 13 therapeutic miRNAs are traced that might safeguard against deregulated lncRNA-mRNA pairs of NORAD, NEAT1, and MIRLET7BHG in FC. This study helps to explain the plausible mechanism of post-COVID neurological manifestation and also to devise therapeutics against it.

Introduction

The whole world is under the mayhem of the coronavirus disease 2019 (COVID‐19) pandemic. Although the lungs are the primary organ for being infected by severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2),[ ] organotropism beyond the respiratory tract, including the kidneys, liver, heart, skin, and brain are also being reported.[ , , ] Viral infection in the brain has been earlier reported to cause havoc. Human coronavirus (HCoV) RNA has been found in the cerebrospinal fluid (CSF) of patients with Parkinson's disease (PD)[ ] and found to be potent in causing neurodegeneration by invading the neurons.[ ] SARS‐CoV‐2 belongs to the family of these coronaviruses and it was also found to be neuro‐invasive.[ ]

In addition to that, the cytokine storm due to SARS‐CoV‐2 infection is massively exhibited by a hyperinflammatory response.[ ] The chronic neuroinflammation induced by cytokine storm is also being treated as a marker of neurodegenerative diseases like Alzheimer's (AD), PD, and Huntington's disease (HD).[ ] During AD progression, pro‐inflammatory cytokines like IL‐1 and IL‐6 inhibit phagocytosis of beta‐amyloid (Aβ) by microglia and so its accumulation causes neuroinflammation.[ ] Whereas, in CSF of PD patients, tumor necrosis factor (TNF) was reported to be responsible for neuronal death in fatal cases.[ ] Intriguingly, neuroinflammation was also found to instigate psychiatric diseases by elicitation of microglial response, altering neuroplasticity, cognition, and behavior.[ ] Thus, aggravation of neuroinflammatory responses by cytokine storm of SARS‐CoV‐2 might cause neurological manifestations in the future.

As per reports from the epicenter of COVID‐19‐ Wuhan, China, the most common neurological manifestation of SARS‐CoV‐2 included dizziness, headache, impaired consciousness, seizures, and acute cerebrovascular disease.[ ] A recent report suggested that a COVID‐19 patient symptomized by fever and abnormal mental status was diagnosed with acute necrotizing encephalopathy which is characterized by blood‐brain barrier (BBB) disruption due to intracranial cytokine storm.[ ] Hence, neurological manifestations are emerging as an aftermath of COVID‐19. Due to the lack of experimental studies in pandemic situation, in silico studies have been priming COVID‐19 research by digging into the utter complexities that might be associated with the comorbidities.[ , ] Several studies with human brain organoids reveal a varying degree of neurotropism for SARS‐CoV‐2.[ , ] Since the brain is a complex heterogeneous tissue consisting of various cell types like–glial, epithelial, and neural cells,[ ] so to understand neurological manifestations of SARS‐CoV‐2 infection, studies should be conducted at the single‐cell level. Yang et al. performed excellent work in studying the dysregulation of coding and non‐coding RNA expression across various cell types in the choroid plexus (CP) and frontal cortex (FC) of patients infected with SARS‐CoV‐2.[ ] However, dysregulation in the transcriptional circuits of these cell types during COVID‐19 has not been elucidated so far. LncRNAs and miRNAs are important classes of non‐coding RNA (ncRNA) that can modulate gene expression post‐transcriptionally. In 2011, Salmena et al. proposed the competing endogenous RNA (ceRNA) hypothesis,[ ] which explained that mRNAs and lncRNAs have one or more common miRNA response elements (MREs). This ceRNA network is being studied for various neurological and psychiatric problems to understand the deregulated regulatory crosstalk between coding and non‐coding RNAs.[ , ] Several types of researches have already been undertaken to unfold the deregulated lncRNA‐mRNA landscape during SARS‐CoV‐2 infection in different organs. It was evidenced that interferon‐mediated inflammatory response is instigated due to anomalous lncRNA‐mRNA interactions during SARS‐CoV‐2 infection of the bronchial epithelium.[ ] Furthermore, cytokine signaling was reported to be aggravated due to deregulated lncRNA‐mRNA interactions in the primary site of COVID‐19, that is, lungs.[ ] Not only lncRNA‐mRNA interactions but also ceRNA networks have already been reported to reveal lncRNAs that might sponge harmful miRNAs in SARS‐CoV‐2 infected lungs.[ ] Thus, to get insights into the malfunctioned regulatory circuits in SARS‐CoV‐2‐infected human brain, we first have reconstructed ceRNA circuits from bulk transcriptomic data of SARS‐CoV‐2 infected CP and FC. Next, we delved into their cell‐specific expression using snRNA‐seq data. Our study infers that downregulation of NORAD, NEAT1, STXBP5‐AS1 in FC and upregulation of MIRLET7BHG and MALAT1 in FC and CP respectively might aggravate serious neurological consequences in the brain. We also have proposed 13 miRNAs which might help to dampen the neurological consequences of COVID‐19.

Results

Reconstruction of ceRNA Networks in SARS‐CoV‐2 Infected Frontal Cortex and Choroid Plexus Tissue of the Brain

We initially identified differentially expressed genes (DEGs) in the FC of SARS‐CoV‐2 infected patients from GSE182297,[ ] where the authors isolated samples from the prefrontal cortex of the brain during autopsies of SARS‐CoV‐2 infected patients. We identified a total of 9254 DEGs by considering a cut‐off of |log2FC|>0.25 and adj. P value < 0.05 (Figure ; File S1, Supporting Information). Then by using human genome annotation in Ensembl Genes 105,[ ] we have obtained 923 DE lncRNAs, 357 DE miRNAs, and 7964 DE mRNAs among the DEGs (Figure 1B). We then isolated 5265 lncRNA–mRNA target pairs among the sets of DEGs by considering their experimentally verified interactions deposited in publicly available datasets‐RISE,[ ] NPinter v 4.0,[ ] and RAID v3.0 (confidence score > 0.5).[ ] With these lncRNA‐mRNA pairs, a network consisting of 3575 nodes and 5265 edges was reconstructed (Figure 1C) in which we have determined 19 hub DE lncRNAs by considering a degree cut off ≥ 10 (Figure 1D). As hub lncRNAs are expected to impact on a greater number of mRNAs in ceRNA networks, so with these hub lncRNAs we have built a ceRNA network. Thus, considering three sets of interactions‐ DE lncRNA‐ DE miRNA, DE miRNA‐ DE mRNA, and DE lncRNA–DE mRNA, we have reconstructed a ceRNA network for FC with 4634 edges and 1100 nodes among which lncRNA, miRNA, and mRNA comprise 12, 84, and 1004 nodes respectively (Figure ).

Figure 1

DE lncRNA‐ DE mRNA network in SARS‐CoV‐2 infected FC of the brain. A) Volcano plot representing the differential expression of genes in FC infected with SARS‐CoV‐2. Genes are considered as significantly expressed only if adj. P value < 0.05 and |log2FC|>0.25. B) Pie chart representing distribution of lncRNAs, miRNAs, and mRNAs among the DEGs in SARS‐CoV‐2 infected FC of the brain. C) DE lncRNA‐ DE mRNA network in SARS‐CoV‐2 infected FC of the brain consisting of 3575 nodes and 5265 edges. D) Table representing 19 hub lncRNAs along with their degrees in DE lncRNA‐ DE mRNA network of SARS‐CoV‐2 infected FC of the brain. Degree cut‐off ≥ 10 was considered for determining hub lncRNAs.

Figure 2

ceRNA network in inSARS‐CoV‐2 infected FC. This network consists of 1100 nodes and 4634 edges among which 12, 84, and 1004 nodes are lncRNA, miRNA, and mRNA respectively.

Similarly, we have identified 3804 DEGs by integrating data from 24 and 72 hpi of SARS‐CoV‐2 in CP organoids from GSE157852[ ] (Figure ,3B; File S1, Supporting Information). Consequently, among these DEGs we have obtained 76 DE lncRNAs, 4 DE miRNAs, and 3732 DE mRNAs (Figure 3C) by using human genome annotation in Ensembl Genes 105.[ ] Here also we have reconstructed a network (Figure 3D) consisting of 523 nodes and 514 edges using 514 DE lncRNA‐ DE mRNA interacting pairs that were identified by following the similar protocol used for FC. From this DE lncRNA‐ DE mRNA network, we have identified 3 DE hub lncRNAs‐ MALAT1, NEAT1, and CASC15 by considering a degree cut‐off ≥ 10 (Figure 3E). Alike FC, we have reconstructed a ceRNA network considering the three sets of interactions. This network consists of 412 edges and 192 nodes comprising of 2 lncRNAs, 3 miRNAs, and 187 mRNAs (Figure ). Thus, these ceRNA networks in FC and CP will help to depict potential cross‐talk between lncRNA ‐miRNA‐mRNA in these regions of the brain.

Figure 3

DE lncRNA‐ DE mRNA network in SARS‐CoV‐2 infected CP of the brain. A) Volcano plot representing the differential expression of genes in 24 hpi of CP organoids with SARS‐CoV‐2 infection. Genes are considered as significantly expressed only if adj. P value < 0.05 and |log2FC|>0.25. B) Volcano plot representing the expression of genes in 72 hpi of CP organoids with SARS‐CoV‐2 infection. Genes are considered as significantly expressed only if adj. P value < 0.05 and |log2FC|>0.25. C) Pie chart representing distribution of lncRNAs, miRNAs, and mRNAs among the DEGs in SARS‐CoV‐2 infected CP of the brain. D) DE lncRNA‐ DE mRNA network inSARS‐CoV‐2 infected CP of the brain consisting of 523 nodes and 514 edges. E) Table representing 3 hub lncRNAs along with their degrees in DE lncRNA‐ DE mRNA network of SARS‐CoV‐2 infected CP of the brain. Degree cut‐off ≥ 10 was considered for determining hub lncRNAs.

Figure 4

ceRNA network in SARS‐CoV‐2 infected CP. This network consists of192 nodes and 412 edges among which 2, 3, and187 nodes are lncRNA, miRNA, and mRNA respectively.

Functional and Pathway Enrichment of the ceRNA in SARS‐CoV‐2 Infected Choroid Plexus and Frontal Cortex Tissue of the Brain

Since downregulation of essential mRNAs or upregulation of harmful mRNAs due to deregulation of lncRNAs involved in the ceRNA network could be fatal, thus we have separately performed functional and pathway enrichment analysis of downstream mRNAs of up and down‐regulated hub lncRNAs in ceRNA axes by integrating results from clusterProfiler v4.0.0[ ] and enrichR v3.0[ ] packages in R v4.1.1. We have considered an adj. P value <0.05 as a cut‐off[ ] to determine significantly enriched functions as well as molecular pathways. Figure Adelineatesthat downregulation of hub lncRNAs in SARS‐CoV‐2 infected FC might downregulate mRNAs enriched for the Biological Process (BP)‐ “cotranslational protein targeting to membrane” [BP: GO:0006613], Molecular Function (MF)‐ “sodium ion transmembrane transporter activity” [MF: GO:0015081] and KEGG Pathway (KP)‐ “GABAergic synapse” [KP: hsa04727] (File S2, Supporting Information). Whereas, upregulation of lncRNAs might upregulate the expression of mRNAs enriched for the BP – “cellular response to decreased oxygen levels” [BP: GO:0036294] and KP – “Proteoglycans in cancer” [KP: hsa05205] (Figure 5A; File S2, Supporting Information). In the case of CP organoids, it was found that upregulation of lncRNAs might downregulate mRNAs involved in the BPs – “cellular response to transform growth factor beta stimulus” [BP: GO:0071560] and “response to nutrient levels” [BP: “GO:0031667”] (Figure 5B; File S2, Supporting Information). These results suggest that deregulation of lncRNAs during SARS‐CoV2 infection hinders the proper expression of mRNAs which may in turn instigate multifarious neurological manifestations in CP and FC tissues of the brain.

Figure 5

Significant function and pathway enrichment of mRNAs that are deregulated in ceRNA networks of A) SARS‐CoV‐2 infected FC and B) SARS‐CoV‐2 infected CP. A cut‐off value for adj. P value < 0.05 is considered and a maximum of five functions/pathways with the lowest adj. P values in each case are represented.

Reconstruction of Cell Specific DE lncRNA‐DE mRNA Network in SARS‐CoV‐2 Infected Choroid Plexus and Frontal Cortex Tissue of the Brain

Our results from bulk RNA‐Seq data have shown that deregulation of hub lncRNA‐mRNA interactions in ceRNA networks of SARS‐CoV‐2 infected brain might incur serious consequences. But, the human brain is a heterogeneous tissue with various types of cells, for example, glial, epithelial, neuronal cells, etc.[ ] To gain detailed insights into the cell‐type specific deregulation of regulatory pathways in virus‐infected CP and FC tissue of the brain, we have analyzed the publicly available snRNA‐Seq dataset GSE159812.[ ] From this, we have identified DE lncRNAs and DE mRNAs from different cell types of CP and FC. However, DE miRNAs were unavailable in snRNA‐Seq data. Using an unsupervised clustering method and based on previously identified markers,[ ] we have classified the CP transcriptomic data into 7 cell types, namely‐ Epithelial, Mesenchymal, Neuronal, Ependymal, Monocyte, Glial, and Endothelial cells. Similarly, we have identified 7 different cell types—Oligodendrocyte, Astrocytes, Oligodendrocyte Progenitor Cells (OPCs), Excitatory Neuron, Inhibitory neuron, Microglia, and Endothelial cells in FC (Figures S1 and S2, Supporting Information). By considering a cut‐off of |log2FC|>0.25 and adj. P value < 0.05, we have identified a total of 392 and 929 DEGs in CP and FC respectively (File S3, Supporting Information). Based on human genome annotation in Ensembl Genes 105,[ ] we have obtained a total of 16 DE lncRNAs and 380 DE mRNAs in different cell types of CP and 60 DE lncRNA along with 869 DE mRNA in different cells of FC (Figure ). With these DE lncRNAs and DE mRNAs we have created cell‐specific DE lncRNA‐DE mRNA networks, as well as identified hub DE lncRNAs, in SARS‐CoV‐2 infected CP and FC by following a similar protocol discussed in the earlier section. After that, to identify the specific cellular location of DE hub lncRNA‐DE mRNA pairs in FC and CP tissues retrieved from bulk RNA‐Seq data, we have mapped them with that of the snRNA‐seq dataset. Mapping of these pairs has revealed only three downregulated lncRNAs—NORAD, NEAT1, and STXBP5‐AS1 which are common in FC and two upregulated lncRNAs‐ MIRLET7BHG and MALAT1 respectively in FC and CP. Subsequently, we have found 7, 4, 3, and 1 inversely correlated NORAD‐mRNA pairs respectively in excitatory neurons, OPCs, Inhibitory neurons, and Oligodendrocytes of the FC (Figure , 7E). In vitro studies have reported that NORAD has neuroprotective roles against PD[ ] and it protects PD substantia nigra from oxidative stress.[ ] Thus, reduced expression of NORAD might cause neurological consequences. Meanwhile, overexpression of NEAT1 was evidenced to be neuroprotective against PD‐mediated oxidative stress[ ] and thus, mutation of NEAT1 showed abnormal reaction to physiological stress.[ ] However, we have noticed that SARS‐CoV‐2 infection causes downregulation of NEAT1 in our datasets. We have found that downregulation of NEAT1 could result in upregulation of SRRM3 in Inhibitory neurons (Figure 7B). Moreover, 4, 2, and 1 positively correlated NEAT1‐mRNA pairs were attained respectively in Astrocytes, Inhibitory neurons, and Oligodendrocytes of FC (Figure 7B, 7D, 7E). For the downregulated STXBP5‐AS1, we have received 10 inversely correlated and 21 positively correlated lncRNA‐mRNA pairs in excitatory neurons (Figure 7A) whereas 2 inversely correlated and 7 positively correlated pairs were observed in inhibitory neurons (Figure 7B). The lncRNA, STXBP5‐AS1 was also reported to be downregulated during tumorigenesis in glioma.[ ] Apart from these 3 downregulated hub lncRNAs, MIRLET7BHG was observed to be upregulated in virus‐infected FC. Moreover, 23 inversely correlated MIRLET7BHG‐mRNA pairs and 8 positively correlated MIRLET7BHG‐mRNA pairs were found in excitatory neurons of FC (Figure 7A). Literary evidence suggested that MALAT1 is involved in inflammasome activation and reactive oxygen species (ROS) generation during the progression of PD.[ ] Thus, the upregulation of MALAT1 observed in CP might be related to the reported phenomenon during SARS‐CoV‐2 infection. Moreover, 8 positively related and 11 inversely related MALAT1‐mRNA pairs were detected in neuronal cells of CP (Figure 7F). Together these results indicate that deregulation of these lncRNA ‐mRNA interactions might be correlated with serious neurological consequences in different cell types of CP and FC.

Figure 6

UpSet plots showing the distribution of DE lncRNAs and DE mRNAs along with the number shared among various cell types of SARS‐CoV‐2 infected CP and FC. A) Distribution of DE lncRNAs in CP. B) Distribution of DE mRNAs in CP. C) Distribution of DE lncRNAs in FC. D) Distribution of DE mRNAs in FC.

Figure 7

Cell specific deregulated lncRNA‐mRNA network due to SARS‐CoV‐2 infection in A) Excitatory Neurons of FC, B) Inhibitory Neurons of FC, C) OPCs of FC, D) Astrocytes of FC, E) Oligodendrocytes of FC, and F) Neuronal cells of CP.

Reconstruction of Probable Cell‐Specific ceRNA Axes in SARS‐CoV‐2 Infected CP and FC of the Brain to Trace Therapeutic miRNAs

Using the snRNA‐Seq dataset (GSE159812), we have received 5 lncRNA mediated subnetworks in different cell types of the brain. However, this dataset is devoid of miRNA expression data. Thus to reconstruct probable cell‐specific ceRNA axes, we have considered significantly DE miRNAs from bulk transcriptomic data, that is, GSE182297[ ] for FC and GSE157852[ ] for CP and mapped them between cell‐specific hub lncRNA‐mRNA pairs. For each downregulated lncRNA in different cell types of FC we have obtained every possible combination of ceRNA axes, that is, down lncRNA‐ down miRNA‐ down mRNA, down lncRNA‐ up miRNA‐ down mRNA, down lncRNA‐ down miRNA‐ up mRNA, down lncRNA‐ up miRNA‐ up mRNA. Cell‐wise ceRNA subnetwork is delineated in Figure . Likewise, for two upregulated lncRNAs we have received the same pattern of combination of mRNA and miRNA with up lncRNA in neuronal cells of FC and CP (Figure 8A,8F).

Figure 8

Cell specific deregulated ceRNA networks due to SARS‐CoV‐2 infection in A) Excitatory Neurons of FC, B) Inhibitory Neurons of FC, C) Astrocytes of FC, D) Oligodendrocytes of FC, E) OPCs of FC, and F) Neuronal cells of CP.

In FC, excitatory neurons seem to be most affected as deregulation of ceRNA axes due to downregulation of NORAD and upregulation of MIRLET7BHG alter the expression of 35 and 11 mRNAs, respectively (Figure 8A). Again, deregulation of ceRNA axes due to downregulation of NEAT1 maximally affects astrocytes by downregulating 4 mRNAs in this cell type. (Figure 8C). In CP, neuronal cells might be at risk because upregulation of MALAT1 might deregulate ceRNA axes that alter the expression of 12 mRNAs (Figure 8F).

Now, it is apparent that miRNAs act as rheostats that fine‐tune gene expression. Thus, downregulated miRNAs in the ceRNA axes where the downstream mRNAs are upregulated could be targeted for designing miRNA‐based therapeutics as an optimal expression of those miRNAs would be able to restore the proper expression of mRNAs. After searching such miRNAs in the ceRNA axes (Figure 8A–F) composing with lncRNA (up/down)‐miRNA (down)‐mRNA(up), we have received13 miRNAs namely‐hsa‐let‐7c‐5p, hsa‐miR‐103b, hsa‐miR‐3139, hsa‐miR‐601, hsa‐miR‐623, hsa‐miR‐125a‐3p, hsa‐miR‐3685, hsa‐miR‐4260, hsa‐miR‐125a‐5p, hsa‐miR‐1284, hsa‐miR‐3909, hsa‐miR‐5191, and hsa‐miR‐4254 from33 ceRNA axes (Figure 8).Further experimental validation is required for these miRNAs to be used as a therapeutic measure against neurological manifestations exerted by SARS‐CoV‐2.

Discussion

The amplification of inflammatory response caused by cytokine storm is the primal response of the body against SARS‐CoV‐2 infection.[ ] But similar forms of immunological chaos have also been found to cause neurodegeneration in multiple sclerosis (MS), PD, AD, and HD.[ ] Neuroinflammation was also reported to be the key factor in psychiatric diseases like bipolar disorder and schizophrenia.[ ] Thus, it was hypothesized that aggravation of neuroinflammatory response might pave the pathway for neurological complications in COVID‐19 patients.[ ]

To proceed with our analysis, we initially screened DEGs from bulk RNA‐Seq datasets (Figure 1A and File S1, Supporting Information) and discovered several important DEGs that have previously been linked to neurological disorders, such as CCL19, which is an upregulated DEG in FC and was found to be expressed in CSF during neuroinflammation in MS,[ ] and IGKV1‐9, encoding the V region of the immunoglobulin light chain that aids in antigen presentation[ ] and thus may be responding to viral infection. In the case of CP (Figure 3A and File S1, Supporting Information), we found upregulated DEGs like CCL7 and IL18 which are important messengers of the immune response as well as key players in neuroinflammatory response.[ , ] Deregulation of these genes indicates that an inflammatory response might be elicited due to SARS‐CoV‐2 infection in both FC and CP. To decipher the dysregulation of non‐coding RNAs associated regulatory circuits, we have reconstructed ceRNA networks in both CP and FC considering three types of interactions‐ lncRNA‐mRNA, miRNA‐mRNA, and lncRNA‐miRNA. For constructing the ceRNA networks, we first reconstructed lncRNA‐mRNA networks and identified hub lncRNAs (Figures 1C and 3D). Then, the ceRNA networks were built focusing on hub lncRNA‐miRNA‐mRNA interactions (Figures 2 and 4). Among the hub lncRNAs in FC, KCNQ1OT1, and SLC25A25‐AS1 are found to be upregulated (Figure 1C). Earlier, KCNQ1OT1 was reported to be overexpressed in traumatic brain injury (TBI) and it stimulated overexpression of cytokines.[ ] Whereas, SLC25A25‐AS1 was evidenced to be highly expressed in males of AD brain.[ ] Interestingly, we have found SYN1 as an interacting mRNA of KCNQ1OT1 which is downregulated in FC. An earlier study had reported that mutated SYN1 facilitates loss of synaptic activity in autism and partial epilepsy.[ ] Similarly, VDAC1, an interacting mRNA of SLC25A25‐AS1, is also found to be downregulated in FC, which resembles a characteristic of the AD brain.[ ] In the case of CP (Figure 3D), the upregulated hub lncRNACASC15 was reported to be involved in lncRNA‐miRNA crosstalk during acute ischemic stroke.[ ] Moreover, one of its downstream mRNACDKN1Awas evidenced to be involved in complement and coagulation cascades in glioblastoma.[ ] Thus, it could be inferred that deregulation of several lncRNA‐mRNA interactions due to viral infection might aggravate neurological disorders in FC and CP.

Several important axes in ceRNA networks that could be implicated for neurological manifestations were also identified (Figures 2 and 4). The upregulated miRNA hsa‐miR‐106b‐5p interlinking the down NORAD‐ down RFX3 in FC (Figure 2) might cause psychiatric problems because hsa‐miR‐106b‐5p was reported to be upregulated in MS patients[ ] and mutated RFX3 was reported in brains with autism spectrum disorder (ASD).[ ] Moreover, upregulated miRNA hsa‐miR‐125b‐1‐3p, which could bind with up MALAT1‐down MBD5 pair in CP (Figure 4), was reported to induce inflammation and oxidative stress in AD models.[ ] Whereas, mutation in MBD5 was found to be associated with cognitive and intellectual disability (ID).[ ] Thus, continued upregulation of this miRNA might cause an inflammatory response in the BBB and thereby lead to mental retardation in the future. Functional and pathway enrichment analysis showed that in FC (Figure 5A), downregulation of hub lncRNAs might downregulate a set of mRNAs that are functionally enriched for sodium ion transmembrane transporter activity. Whereas, upregulation of hub lncRNAs might upregulate the expression of mRNAs that are enriched for the function‐ cellular response to decreased oxygen levels. These results indicate SARS‐CoV‐2 infection might be instigating hypoxia‐mediated neuronal death by causing impaired ion channel response in FC.[ ] In the case of CP (Figure 5B), upregulation of hub lncRNAs might downregulate mRNAs that are functionally enriched for nutrient level sensing. Thus, deregulation of these mRNAs could deprive the CP of adequate nutrients required for its proper functioning. The results of bulk RNA‐Seq analyses stirred us to find out the cellular location of deregulated pairs participating in regulatory networks in FC and CP tissues of the brain.

Subsequently, mapping of hub lncRNA‐mRNA pairs from bulk RNA‐Seq data with that of snRNA‐Seq data (Figure 7) and reconstruction of cell‐specific ceRNA networks by tracing miRNAs from bulk RNA‐Seq that interlink cell‐specific hub lncRNA‐mRNA pairs showed us that two hub lncRNAs NORAD and NEAT1 are downregulated and one hub lncRNA MIRLET7BHG is upregulated in various cell types of FC (Figures 7 and 8). Literary evidence suggested that a functional NORAD could prevent apoptosis and mitochondrial dysfunction in PD.[ ] Although, the absence of NORAD was not evidenced to directly elicit an inflammatory response but the improper clearing of apoptotic cells by phagocytes could result in hyper‐inflammation.[ ] Moreover, we have found that NORAD is positively correlated with the mRNA TOMM20 in excitatory neurons of FC (Figure 7A). It was reported that TOMM20 encodes a non‐functional mitochondrial import protein which play role in pathogenesis of PD.[ ] At the same time, we have noticed that the downregulation of NORAD might downregulate PEG10 and we traced 6 downregulated miRNAs interlinking this axis (Figure 8A). Experiments in mouse models demonstrated that PEG10 might regulate mood, emotion, and circadian rhythm in the brain.[ ] Thus, consistent with the earlier reports of neurodegenerative diseases,[ ] we found that downregulated NORAD deregulates the expression of several mRNAs which might have implications for neurological consequences during SARS‐CoV‐2 infection.

NEAT1 was evidenced to be both anti‐apoptotic and anti‐inflammatory in case of traumatic brain injury.[ ] Here, we have noticed downregulation of NEAT1 might be causing upregulation of TBL1XR1 in astrocytes (Figure 7D). It was evidenced that mutation of TBL1XR1 causes behavioral abnormalities with delayed motor functions in humans.[ ] Again downregulation of NEAT1 might downregulate HIVEP2 in inhibitory neurons and astrocytes (Figure 7B,7C). Interestingly, loss of HIVEP2 was reported to cause ID.[ ] This cell‐specific variability in expression shows that proper neural transmission and normal behavioral anomalies might be exacerbated due to the downregulation of NEAT1 in SARS‐CoV‐2 infected brains. Upregulation of MIRLET7BHG was found to be inversely correlated with mRNAs ARPP19 and DYNC1H1 in excitatory neurons (Figures 7A and 8A). Low levels of ARPP19 were reported to be the characteristics of the cortical region in AD brain.[ ] Also, mutations in DYNC1H1 encoding a heavy chain of dynein motor protein, are associated with several neurological diseases.[ ] Thus, upregulation of MIRLET7BHG might be deregulating several mRNAs which might culminate in aggravation of serious neurological consequences.

In the case of CP (Figure 7F), we found that MALAT1 is upregulated in both bulk and neuronal cells of the snRNA‐Seq dataset. Downregulation of MALAT1 was evidenced to reduce neuronal cell death[ ] and it was also reported that MALAT1 could promote neuroinflammation in PD mouse models.[ ] Thus, upregulation of MALAT1 in neuronal cells could be fatal. We found that upregulation of MALAT1 might downregulate the expression of SOX5. Loss of function of SOX5 was found to be associated with amyotrophic lateral sclerosis.[ ] Again SEMA6D was also found to be downregulated due to the upregulation of MALAT1. SEMA6D was reported to be associated with synaptic plasticity and the loss of this gene might affect reading ability in children.[ ] Thus, upregulation of MALAT1 might also cause neurological disorders in patients infected with SARS‐CoV‐2.

In the end, we aimed to trace miRNAs that could be used therapeutically to target deregulated molecular pathways as miRNA‐based therapeutics are evidenced to be promising.[ ] It could be speculated that in the ceRNA axes, where deregulated lncRNA–downregulated miRNA‐ upregulated mRNA is observed to occur, there might be a chance to restore optimal expression of these mRNAs by increasing expression of the miRNAs. Keeping this in mind, we have focused on these types of ceRNA axes and found that 13 miRNAs were downregulated along 33 deregulated ceRNA axes in which downstream mRNAs are upregulated (Figure 8). Interestingly, some of these upregulated mRNAs had already been reported to be correlated with neurological diseases like PLD1, SLC7A2, and SH3PXD2A. High expression of PLD1 and SLC7A2 was reported to be a characteristic of AD brain[ , ] while SLC7A2 was also found to be upregulated in human HD brain.[ ] Moreover, an earlier study had depicted that SH3PXD2A is associated with ADAM12 to enhance neurotoxicity during the progression of AD.[ ] Whereas, most of these 13 downregulated miRNAs were already reported to be useful for the brain, for example, hsa‐mir‐125a was found to interfere with the viral translation process.[ ] On the other hand, hsa‐let‐7c‐5p was reported to play a neuroprotective role by preventing microglial activation during cerebral ischemia injury.[ ] As microglia activation might trigger cytokine storm in COVID‐19 infection,[ ] thus this miRNA might be therapeutically important against SARS‐CoV‐2 infection of the brain. hsa‐miR‐103b has already been proposed as a therapeutic miRNA that could promote neurite outgrowth in AD models.[ ] Moreover, hsa‐miR‐124‐3p was evidenced to reduce inflammatory response by inhibiting activation of microglia and astrocytes.[ ] Thus, we conclude that optimal expression of these 13 miRNAs might be a good therapeutic measure to safeguard against 33 deregulated ceRNA axes in SARS‐CoV‐2 infected brain. Further experiments need to be carried out in the future to unravel the efficacy of these miRNAs in combating the disease. Moreover, the study is limited because in case of bulk RNA‐seq dataset for FC, the authors from the original study[ ] considered two technical replicates from the prefrontal cortex of a SARS‐CoV‐2 infected patient. So, our findings need to be validated by integrating similar bulk RNA‐seq datasets when they are generated in future, as more than one biological replicate could reduce sample‐specific biasness in the data.

Conclusion

In this exploratory study, we have represented a landscape of deregulated lncRNA‐miRNA‐mRNA interaction which could be correlated with neurological complications during or after SARS‐CoV‐2 infection. Analyzing Bulk and snRNA‐Seq data, we have deciphered the roles of four hub lncRNAs‐ NORAD, NEAT1, STXBP5‐AS1, and MIRLET7BHG in aggravating neurological consequences via deregulating several downstream mRNAs in different cell types of FC. Moreover, in CP, MIRLET7BHG, and MALAT1 were found to be crucial for the same purpose. We have also depicted 13 miRNAs in FC which could target 13 over‐expressed harmful mRNAs in 33 ceRNA axes. Targeted expression of these miRNAs might help to cope with the adverse condition in the brain during infection. Moreover, the expression of miRNAs was not available in the snRNA‐seq dataset. So, cell‐specific experiments should be conducted to elucidate their therapeutic role in the future. The study is also limited in terms of RNA‐seq datasets used. As only one bulk RNA‐seq dataset for each condition was available during the study, integrating similar datasets in the future would further validate our findings. This in silico study will help in understanding cell‐specific variability in ncRNA interactions and design therapeutic measures accordingly.

Experimental Section

Data Processing and Quality Check

Initially, the raw reads were downloaded for bulk RNA‐seq data from GSE182297,[ ] GSE157852,[ ] and snRNA‐Seq data in GSE159812[ ] available in NCBI (National Center for Biotechnology Information). GSE182297[ ] contained 2 brain controls and 2 samples from prefrontal cortex of patients infected with SARS‐CoV‐2. GSE157852[ ] contained 3 control samples and 3 test samples each for 24 hours’ post‐infection (hpi) and 72 hpi of SARS‐CoV‐2 in CP organoids. Basic quality checks with the raw counts were done to remove outliers by visualizing samples in R v4.1.1. For snRNA‐Seq data GSE159812,[ ] 6 controls and 4 SARS‐CoV‐2 infected samples were used from lateral ventricles of CP. Similarly, 4 control and infected samples were considered from parenchymal cells of FC. The sample‐specific patient details are mentioned in File S4, Supporting Information. The read counts were processed using Seurat v4.0[ ] in R v4.1.1. After visually examining the counts, cells were filtered that have unique feature counts >200 or <4000 and <25% mitochondrial counts for CP. Similarly, for FC, cells having unique feature counts >200 or <5000 and having <5% mitochondrial counts were considered. Next, principal component analysis (PCA) was performed on the normalized data after limiting highly variable features to 2000. After PCA, the samples were integrated using RunHarmony[ ] and then used 30 dimensions as input for Seurat's RunUMAP, FindNeighbors, and FindClusters (at 0.2 resolution) functions for both the tissues. The Louvain algorithm score for cluster identification was 0.92 and 0.97 for CP and FC respectively. The cell clusters were identified by considering positive differential expression of previously devised markers[ ] of each cluster against all other clusters using MAST algorithm[ ] (Figures S1 and S2, Supporting Information).

Identification of Differentially Expressed Genes from Bulk and Single‐Cell Transcriptome Data

For identifying DEGs in SARS‐CoV‐2 infected FC, bulk transcriptomic data from GSE182297 were used, where the authors collected two technical replicates from the prefrontal cortex of a SARS‐CoV‐2 infected patient (File S4, Supporting Information). For CP, transcriptomic data from SARS‐CoV‐2 infected CP brain organoids in GSE157852[ ] was used. For this dataset, an exhaustive list of DEGs was prepared by combining DEGs from both 24 and 72 hpi. The read counts were analyzed using DESeq2[ ] in R v4.1.1. A cutoff of |log2FC| > 0.25 and adj. P value <0.05 was considered for identifying DEGs.

Once the cell type specific clusters were identified in snRNA‐Seq data, DE lncRNA and DE mRNAs from different cell types of CP and FC were then identified. For this MAST[ ] was used to compute differential expression within a specific cell type from SARS‐CoV‐2 infected CP and FC. To reduce sample‐specific biasness during DEG identification, a generalized linear mixed model using zlm() function was fitted with a random effect for samples in case of both the tissues. The controls and COVID‐19 samples for each cluster were compared for cluster‐specific DE RNA identification. A cut‐off of |log2FC|≥0.25 and adj. P value <0.05 were used to screen DEGs. Pearson correlation coefficient (PCC) analysis was also performed to confirm that the identification of cell‐specific DEGs was not biased by number of nuclei isolated for each cell type (File S3, Supporting Information). Then human genome annotation (GRCh38.p13) available in Ensembl Genes 105[ ] was used to segregate DEGs into DE lncRNA and DE mRNA.

Reconstruction of Cell Specific Hub DE lncRNA‐ DE mRNA Pairs, and ceRNA Networks

For identifying cell‐specific DE lncRNA‐ DE mRNA pairs, common interactions were used from three lncRNA‐ mRNA interaction databases‐ RISE,[ ] NPinter v 4.0,[ ] and RAID v3.0 (confidence score>0.5).[ ] These common sets of interactions were uploaded to Cytoscape 3.8.2. after which network analyzer was used to determine hub DE lncRNAs. To construct ceRNA network, miRNAs that interact with both hub DE lncRNA and its partner DE mRNA were identified by considering lncRNA‐miRNA interactions from RAID v3.0(confidence score>0.5)[ ] and miRNA‐mRNA interactions from miRWalk v3 (binding probability = 1).[ ] Both 5'UTR and 3'UTR miRNA‐mRNA interactions were considered in the study. All the interaction networks were visualized in Cytoscape 3.8.2.[ ]

Functional and Pathway Analysis of Cell Specific lncRNA‐mRNA Networks

For GO (Gene Ontology) Biological Process (BP), GO MF, and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis, clusterProfiler v4.0.0[ ] and enrichR v3.0[ ] packages were used in R v4.1.1 by considering an adj. P value <0.05. For GO functional and pathway enrichment analysis, hypergeometric test (one‐tailed variant of Fisher's exact test) was used to identify significant over‐represented GO terms and pathways.

Statistical Analysis

For GSE182297,[ ] 2 brain controls and 2 samples from prefrontal cortex of patients infected with SARS‐CoV‐2 were considered. For GSE157852[ ] 3 control samples and 3 samples each for 24 and 72 hpi of SARS‐CoV‐2 infection in CP organoids were considered. For snRNA‐seq data GSE159812,[ ] 6 controls and 4 COVID‐19 infected samples were used from CP and 4 control and infected samples from parenchymal cells of FC. While processing the snRNA‐Seq data, PCA was performed and considered 30 PCs for linear dimensional reduction on the data, UMAP (Uniform Manifold Approximation and Projection) for non‐linear dimensional reduction, and Louvain algorithm for cluster identification. The DEGs were identified by determining fold change cut off of |log2FC|>0.25 and adj. P value < 0.05 across all RNA‐Seq and snRNA‐Seq data. For GO functional enrichment analysis and pathway enrichment analysis hypergeometric test, which was one‐tailed variant of Fisher's exact test, was used to identify significant over‐represented GO terms and pathways. An adj. P value < 0.05 was considered as cut‐off for functional and pathway enrichment analysis.

Conflict of Interest

The authors declare no conflict of interest.

Supporting Information

Click here for additional data file.

Supporting Information

Click here for additional data file.