Astrocyte infection, RNA-seq, bioinformatic analysis, and validation of results

Human astrocytes (primary cells) were infected at a low passage number (4th passage) to prevent any genetic and phenotypic changes in the cells that could affect gene expression. It should be noted that the cells were obtained commercially, and it was not possible to determine whether the astrocytes were of male or female origin. Astrocytes were challenged with the SARS-CoV-2 variants with the multiplicity of the infection 0.2. The low MOI was preferred in our investigation since several studies have reported low viral load in the brain and cerebrospinal fluid of COVID-19 patients^12,13,14,15, indicating a low number of the virions that potentially infect the brain cells during the early phase of neuroinfection. The low MOI was also preferred considering the restricted number of pathogens that can get across the impermeable blood-brain barrier and glia limitans, i.e., from the blood to the brain parenchyma^15,16,17. Please note that various MOIs and incubation time-points are found in the literature. For instance, astrocytes and brain microvascular endothelial cells were infected with a low MOI of 0.01 to 0.1 for 6 to 24 h in several studies^18,19,20,21, the BBB model comprised of endothelial cells was infected with an MOI 1 and 10²², whereas, the brain organoids were infected with an MOI of 0.01 for three days²³. Although the low MOI employed in our study may aid in studying early infection events in a controlled manner, it should be noted that it may not provide a comprehensive picture of cell signalling events following viral replication in late stages of infection, which call for spatial transcriptomic analysis using varying MOIs.

Nine cDNA libraries were generated from three biological replicates: non-induced astrocytes (libraries NK1 to NK3), cells induced with strain Slovakia/SK-BMC-BA42/2022 Omicron variant (A1 to A3), and strain Slovakia/SK-BMC-BA15/2021 Delta variant (B1 to B3) (Supplementary information 1, Table S1). The thorough quality control of data obtained from NGS showed the sequencing depth varies from ~ 9.7 million (NK2) to ~ 10.62 million (B2) reads. Complete details of quality control, alignment and splices, results of the mapped regions, the strand specificity of the sequenced reads, read assignment to genes (single-end reads were aligned to the human reference genome), biotypes of the assigned genes, differential gene expression analysis, and its normalization are presented in supplementary information 1. The percentage of uniquely and multi-mapped reads was high for all samples (~ 77% of uniquely mapped and ~ 12% of multi-mapped), and the total number of mapped reads varied between ~ 9.6 M and ~ 11 M (Supplementary Information 1 Table S2, Figures S1 and S2). The mapping distribution (normalized gene body coverage) and library strand specificity are presented in supplementary figures S3 and S4. Total number of mapped genes are in supplementary information 2 (dataset 1). RNA-seq analysis was validated by qRT-PCR with a subset of representative genes. Correlation between RNA-seq and qRT-PCR was r = 0.87 (P 1 panel A), indicates reliability of RNA-seq analysis data.

Differentially expressed genes (DEGs)

The maximum and minimum log2fold change values (logFC), resulting from averaging three independent biological replicates, are presented in Fig. 1 panel B. In this study, a gene was classified as differentially expressed (DEG) if its expression exceeded the cutoff +/- 0.9. Analysis of RNA-seq expression profiles performed with edgeR of Bioconductor revealed a total of 346 DEGs (197 up-regulated genes and 149 down-regulated) evoked in astrocytes challenged with the Omicron variant, while 341 DEGs (215 up-regulated genes and 126 down-regulated) were found in cells infected with Delta (Supplementary information 2, datasets 2 and 3). A total of 215 genes were evoked in a similar fashion in both treatments (149 up-regulated, 66 down-regulated), whereas one gene, SNORD14E, affiliated with the snoRNA class, was upregulated in astrocytes infected with Omicron and downregulated by the Delta variant (Supplementary information 2, dataset 3). A surprisingly significant number of the genes were exclusively evoked by Delta (82 upregulated and 48 downregulated) and Omicron variants (65 upregulated and 60 downregulated, Fig. 1 panel C). The principal component analysis (PCA) performed on the DEGs revealed consistent clustering between biological replicates for each condition (Fig. 1D, Supplementary information Figure S5). Thus, it can be concluded that the differences are emerging from the treatment and not due to the discrepancies between the replicates.

It is important to note that in transcriptome analysis, the level of gene expression and pattern of expression may vary depending on the length of incubation. Six hours of the incubation was performed in the present study to obtain comprehensive data of the genes evoked in the early phase of infection.

Categorization of the DEGs according to the gene ontology (GO)

Differentially expressed genes were categorized based on their involvement in different biological processes using a peer-reviewed Reactome server that uses enrichment analysis corrected for false discovery rate using the Benjamani-Hochberg method (https://reactome.org). A complete list of the GO-biological processes and pathways for each treatment found in the enrichment analysis is presented in supplementary information 2 datasets 4 and 5, while the top 20 up- and downregulated genes are listed in datasets 6 and 7. Pathway enrichment analysis using the Reactome server revealed significantly altered expression of five genes in infected astrocytes: BST2, TLR2, PARP16, SERPINE1, and CHMP4C (Fig. 1 panel E), all involved in the viral infection pathway (SARS-CoV-2 infection pathway).

Top 20 DEGs (DESeq2, with independent filtering, adj.p-value = +/-0.9) when compared between treatments and among replicates showed clear distinction between the gene expression in the control cells and challenged with SARS-CoV-2 variants (Fig. 2 panel A), indicating that the alterations in the gene expression (DEGs) are not due to variations among the replicates. The pathway enrichment analysis performed with g: Profiler (g: GOst, with data source Reactome pathways) enriched DEGs into pathways, and the top ten pathways activated by each virus were signal transduction, neuronal system, transcription, signalling via GPCR, transcription, immune system, etc. (Fig. 2, panel B).

Signaling events evoked in astrocytes

The virus attachment on the host cell surface emanates a series of events in the cells of the neurovascular unit, namely viral sensing by pattern-recognizing receptors, induced virus uptake, and transcytosis. Attachment can also evoke induced reorganization of the extracellular matrix and dysregulation of cell-cell junction molecules (adherens, gap, and tight junctions, Fig. 2 Panel C). Initial stages of the viral infection are often associated with disturbance in the cell physiology (increased level of stress proteins), altered expression of cell-cell signaling molecules, and the activation of the innate immune system. In the case of astrocytes, altered expression of neurotransmitters is also of particular importance. Below we have elaborated on biological events that we observed in the astrocytes in response to Delta and Omicron variants and compared our results with available literature. We touch on five major biological processes that could be related to neurological disorders often seen in COVID-19 patients. These processes are (1) PAMP (pathogen-associated molecular pattern) recognition by astrocytes and cytokine response, (2) neurotransmitter signaling, which may affect astrocytic and neuronal functions, (3) dysregulation in the neurotrophins, (4) changes in the expression of proteins that form various cell-cell junctions in the blood-brain barrier and glial barrier (glia limitans; Fig. 2 Panel C), and (5) change in non-coding RNA expression.

PAMP recognition by astrocytes and cytokine response

SARS-CoV-2 sensing by pattern recognition receptors

The innate immune system can recognize pathogens via pattern recognition receptors (PRRs). Six types of PRRs can sense the molecular patterns of pathogens, which are the toll-like receptors (TLRs), C-type lectin receptors (CLRs), nucleotide-binding and oligomerization domain (NOD)-like receptors (NLRs), absent in melanoma 2 (AIM2)-like receptors (ALRs), retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs), and cyclic GMP-AMP (cGAMP) synthase (cGAS). Several PRRs have been reported to be involved in sensing β-coronavirus infection, including melanoma differentiation-associated protein 5 (MDA-5, IFIH1)²⁴, NLR family pyrin domain containing 3 (NLRP3, NLRP3), TLR7 (TLR7) and TLR2 (TLR2)^25,26. Thus, it could be predicted that among those PRRs, several receptors could also recognize SARS-CoV-2 infection, causing initiation of an intracellular signaling cascade. So far, activation of the TLRs, RLRs, NLRs, ALRs, and cGAS, signaling pathways in response to SARS-CoV-2 were documented in different cell types²⁷. In particular, the TLR2 signaling pathway was induced by E protein of SARS-CoV-2²⁸ and TLR4 was activated in response to the S protein²⁹. In our data, the astrocytes induced by Delta and Omicron variants, we could find significant dysregulation of TLR2 by Omicron (logFC 0.99; Fig. 2, panel D), while in case of Delta, the upregulation was seen only by 0.65 log fold. None of the other TLRs were significantly activated.

The pathway enrichment analysis of our data revealed that 5 genes involved in C-type lectin receptor pathway were evoked (Fig. 2, panel D), in which the Fc Epsilon Receptor Ig encoded by FCER1G was substantially downregulated in both Delta (logFC − 1.19) and Omicron (logFC − 2.42) infected astrocytes. Expression of the other members of the CLR family, such as L-SIGN (CLEC4M) and LSECtin (CLEC4G), remained unaltered in astrocytes. It is important to note that CLRs are expressed mainly by antigen-presenting cells such as dendritic cells and macrophages, and it was shown that several CLRs (CLEC10A, ASGR1, LSECtin, L-SIGN, and DC-SIGN) can detect protein S in a glycan-dependent manner³⁰. The recognition of S protein by CLRs causes the production of pro-inflammatory cytokines such as IL-6, which is associated with a cytokine storm in severe cases of COVID-19 infection³⁰. In our dataset transcripts for DC-SIGN and CLEC10A were absent, while expression of ASGR1 was downregulated only by -0.3 log fold in astrocytes infected either by the Delta or Omicron variant (supplementary information 2 dataset 1).

PRRs like NLRs, RLRs and RIG-I were also reported to be involved in the recognition of SARS-CoV-2 molecules. For instance, results from other reports have shown that NLRP3 recognizes GU-rich single-stranded RNA to induce the secretion of pro-inflammatory cytokines in macrophages and is also able to detect S protein, which leads to its induction^31,32. In contrast, in our study, we found slight downregulation of NLRP3 (logFC − 0.11 both variants) in infected astrocytes. RLRs and RIG-I are the key cytosolic sensors of virus RNA. Stimulation of both PRRs have been reported by SARS-CoV-2 in dendritic cells, macrophages or the airway epithelial cells leading to high expression of MDA-5 and RIG-I (DDX58)³³. In astrocytes, however, we found that expression of RIG-I was insignificantly increased (logFC 0.24 Delta, logFC 0.53 Omicron), while MDA-5 expression was slightly downregulated (logFC − 0.2 Delta, logFC − 0.3 Omicron; supplementary information 2 dataset 1). The PRR expression levels found in our investigation, especially for NLRP3, RIG-I, and MDA-5, contradict the expression levels obtained in other cells (e.g., macrophages or microglia), where significant upregulation was reported after viral infection^34,35,36. Although these PRRs are constitutively expressed in astrocytes, the basal expression levels and their stimulation in response to a pathogen assault may differ in astrocytes depending on the region of the brain, animal species, age, etc. (reviewed in³⁷). It is important to note that astrocytes, major cells that maintain homeostasis in the CNS, may exert pro-inflammatory roles during various assaults such as viral infection, but can also exert anti-inflammatory roles to maintain homeostasis and limit inflammation by creating a barrier between lesions and healthy tissue (barrier, also referred to as adventitial cuffs)^37,38. In the CNS, adventitial cuffing is well documented in infections and inflammation, and astrocytes have long been identified to restrict the spread of inflammation away from adventitial cuffs (e.g., into the neural parenchyma) by attenuating the inflammatory pathways³⁹. Earlier intriguing work has showed that only microglia in the brain can activate NLRP3, form a functional NLRP3 inflammasome, and release IL-1β, while astrocytes do not produce IL-1β or induce expression of the NLRP3 inflammasome components in response to similar assaults³⁴. Downregulation of NLRP3 inflammasome associated proteins such as NLRP3 (logFC − 0.11 in both variants), IL-1β (IL1B; logFC − 1.89 Delta, logFC − 1.31 Omicron) and caspase 1 (CASP1; logFC − 0.34 Delta, logFC − 0.24 Omicron) observed in our investigation warrant further research to determine whether astrocytes actually employ different PRRs to sense PAMPs derived from SARS-CoV-2 to trigger downstream pathways.

Cytokine signaling

Following the recognition of PAMPs such as protein S, protein E, nucleic acid, etc., profound dysregulation in gene expression related to cytokine pathways (the chemokine and interferon pathways) is recorded in various cells by a plethora of studies^40,41,42. However, there is limited information available on cytokine signaling events in the cells of the neurovascular unit, particularly astrocytes. Data available so far from a sound experiment reveal significantly altered inflammatory gene expression, reactivity characteristics, and increased cytokine signaling in cortical astrocytes infected with the Wuhan variant⁴³. In contrast to said observations, our results present comparatively low inflammatory and cytokine responses to the Omicron or Delta variant. For instance, only 16 protein-coding genes were evoked in the cytokine signaling pathway (Fig. 2, panel E), in which 9 genes belonged to interferon signaling and 7 genes were from signaling by interleukins (mainly IL-4/IL-13 signaling and IL-10 signaling). The difference in the cytokine response that we observed and the robust response reported by Andrews and colleagues⁴³, could be attributed to different variants used to challenge the cells (i.e., Wuhan used by Andrews and colleagues vs. the Omicron and Delta variants used by us).

It is also important to note that inflammatory reactions could be both beneficial and detrimental to the brain, depending on the strengths of their activation. Mild activation of astrocytes usually reveals neuroprotective effects; however, strong activation gives rise to cytokine dysregulation, which accelerates cognitive impairment and even neurodegeneration⁴⁴. Thus, as we expected, the profound cytokine response in astrocytes was not observed in this study.

Gliotransmitter signaling

Astrocytes respond to physical assault and infections also by releasing transmitters (known as gliotransmitters) that have feedback actions on neurons and actively produce neuroactive molecules such as glutamate, gamma-aminobutyric acid (GABA), D-serine, adenosine 5’ triphosphate (ATP), etc. Since it was recently demonstrated that astrocytes are a suitable site for the multiplication of SARS-CoV-2, disruption of astrocytic functionality may lead to a decrease in neuronal viability, brain disorders, or abnormal behaviour⁴⁵. Changes in glutamate and GABA signaling in astrocytes may be associated with these symptoms⁶. As a result, we sought to investigate how these signaling pathways (glutamate, GABA and purinergic signaling) are altered in astrocytes in response to Delta and Omicron variants.

Glutamate signaling

Neurotransmission of glutamate is mediated by ionotropic (iGluRs) and metabotropic (mGluRs) receptors. We were particularly interested in changes in the expression of major ionotropic receptors, including α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPA), kainic acid receptors (Kainate), N-methyl-d-aspartate receptors (NMDA), and specific glutamate delta-1 (GluD1) and glutamate delta-2 (GluD2) receptors. Among the metabotropic receptors, mGluRs 1–8 were intriguing. Transcriptome analysis revealed disruptions in the expression of several of these receptors (Fig. 3). Among metabotropic receptors, expression of mGluR5 (GRM5) was evoked significantly by both viral variants (logFC 1.86 Delta, logFC 1.12 Omicron). An imbalance in the expression of this receptor is associated with mood swings, anxiety, and learning and memory problems^46,47. It is also hypothesised that mGluR5 upregulation may cause suicidal ideation or mood disorders⁴⁸. mGluR4 (GRM4; logFC − 0.67 Delta, logFC − 0.68 Omicron) and 8 (GRM8; logFC − 0.11 both variants) were slightly downregulated in challenged astrocytes (Fig. 3), while the expression of other mGluRs (1, 2, 3, 6, 7) remained unchanged.

Among ionotropic receptors, the majority of the genes encoding AMPA (GluA1-4, GRIA1-4) were either unchanged or slightly downregulated (such as GRIA1, logFC − 0.16 Delta, logFC 0.01 Omicron; GRIA4, logFC 0.06 Delta, logFC − 0.02 Omicron; Fig. 3); the maximum downregulation observed for GRIA3 logFC was − 0.21 in Omicron-infected cells. It is of note that the gene encoding kainate receptors (GluK1-5) had altered expression patterns, particularly for GluK1 (GRIK1), which was significantly downregulated (logFC − 1.15) following treatment with the Delta variant, whereas GluK2 (GRIK2) was upregulated (logFC 0.84) following the infection with Omicron. Transcriptomic data revealed a greater imbalance in NMDA receptor-related transcripts after infection. Infected astrocytes, for example, showed an increase in GluN2A (GRIN2A) gene expression (logFC 1.3 both variants), while GluN2C (GRIN2C) expression was significantly decreased (logFC − 0.9 for both variants; Fig. 3). Similar to our findings, a previous study reported an imbalance in glutamate receptor expression following infection with neurotropic viruses such as West Nile virus (WNV) and Japanese encephalitis virus (JEV), which resulted in a significant downregulation of mGluR1–5 receptors and NMDA receptors (GluN1 and GluN3A)⁴⁹.

Astrocytes possess the ability to produce glutamate by hydrolyzing glutamine with the glutaminase (GLS) enzyme. The optimal levels of both GLS isoforms, GLS1 (GLS) and GLS2 (GLS2), are critical for brain homeostasis. Upregulation of GLS has been recorded in various brain disorders^50,51. Our data showed an insignificant increase in the expression of GLS1 caused by both variants (logFC 0.2), while expression of GLS2 was significantly downregulated in Omicron variant-infected astrocytes (logFC − 2.03) and remained at basal levels in Delta variant-infected cells (logFC − 0.03; supplementary information 2 dataset 1). Such a pronounced downregulation of GLS was also recorded in other neuroinfections caused by WNV and JEV⁴⁹.

To regulate glutamate levels, astrocytes can either uptake it from the extracellular space or release it through exocytosis. Excitatory amino acid transporters (EAAT) encoded by SLC1A genes regulate uptake, whereas vesicular glutamate transporters (VGLUT) encoded by SLC17A genes regulate exocytosis (Ca2+-dependent or small intracellular vesicle-activated exocytosis)⁵². Our transcriptomic data revealed a significant imbalance in the expression of excitatory amino acid transporters (Fig. 3). Both Delta and Omicron infections significantly reduced the levels of SLC1A3 (logFC − 0.71 Delta, logFC − 1.26 Omicron) and SLC1A2 (logFC − 1.16 Delta, logFC − 1.15 Omicron), which encode EAAT1 and 2, respectively (Fig. 3). SLC1A1 expression, which encodes the EAAT3 transporter, diminished in infected astrocytes, but only marginally (logFC − 0.18 Delta, logFC − 0.55 Omicron), whereas expression of other transporters, such as EAAT4 and 5, was unaffected. Other studies have also reported a decrease in EAAT levels, specifically astrocytic EAAT2, during neuroinfections caused by herpes virus, human immunodeficiency virus 1 (HIV-1), or WNV^53,54,55. Our data also revealed that neither Delta nor Omicron causes imbalance in the expression of vesicular glutamate transporters, except for VGLUT1, whose expression was slightly upregulated in the Omicron-infected astrocytes (logFC 0.67; Fig. 3). In contrast to our findings, a study of mouse brains infected with WNV and JEV found significant downregulation of VGLUT (SLC17A)⁴⁹.

GABA signaling

Astrocytes are capable of producing the GABA and uptaking it from the synaptic cleft⁵⁶. To convert glutamate into GABA, hippocampal astrocytes utilize two isoforms of glutamate decarboxylase GAD67 (GAD1) and GAD65 (GAD2)⁵⁷. In our transcriptomic data, GAD67 had altered expression (logFC 0.39 Delta, logFC 0.75 Omicron), whereas GAD65 remained unchanged (Fig. 4 panel A).

Astrocytes can also produce GABA by oxidizing putrescine with monoamine oxidase B (MAO-B, MAOB). This GABA synthesis pathway involves two additional enzymes, aldehyde dehydrogenase 1 family A1 (ALDH1A1) and sirtuin 2 (SIRT2)⁵⁸. Expression levels of all these enzymes were lowered in Delta and Omicron-infected astrocytes, with a particularly significant downregulation of ALDH1A1 (logFC − 1.0, Fig. 4 panel A) observed after infection with the Delta variant. A recent study, however, reported increased levels of MAO-B in the brains of COVID-19 patients with depressive and reduced cognitive behaviour⁵⁹.

Putrescine can also be converted to GABA through a diamine oxidase (DAO)-dependent pathway, with ALDH1A1 serving as a key downstream enzyme and bestrophin 1 (BEST1) and leucine-rich repeat-containing protein 8 A (LRRC8A) facilitating GABA release⁶⁰. The gene encoding DAO, AOC1, was not found in our RNA-seq data. In general, levels of bestrophins (BEST1-4) were elevated in infected astrocytes, except BEST2, with the highest upregulation observed for BEST4 (logFC 1.0 both variants). The expression of LRRC8A was not significantly evoked by infection with either variant (logFC 0.23 Delta, logFC 0.25 Omicron; Fig. 4 panel A).

Astrocytes also function as GABAceptive cells, which means they express GABA receptors and multiple GABA transporter proteins (GATs)⁶¹. GABA uptake requires the presence of both metabotropic (GABA_AR) and ionotropic (GABA_BR) receptors. The majority of the genes encoding ionotropic receptors changed expression in our study (Fig. 4 panel A), with the most significant decrease observed for GABA_AR beta-2 subunit (GABRB2, logFC − 0.94 Delta, logFC − 1.95 Omicron). This subunit is highly expressed in the brain, and altered expression has been linked to neuropsychiatric or neurological disorders^62,63. The level of another subunit, GABA_AR gamma-3, was significantly increased in our study (GABRG3, logFC 2.58 Delta, logFC 1.18 Omicron; Fig. 4 panel A).

Among GABA transporters, GAT2 (SLC6A13) was the most downregulated candidate, in astrocytes infected with the Omicron variant (logFC − 1.44). Expression levels of other transporters remained near basal level (Fig. 4 panel A). Numerous neurological illnesses, such as depression and anxiety⁶⁴, autism⁶⁵, schizophrenia, and bipolar disorder⁶⁶, have already been linked to defects in GABA signaling. This highlight needs more research in this area to determine the relationship between GABA signaling pathway disorder and neurological disorders seen in COVID-19.

Purinergic signaling

In response to viral infection, a variety of cell types, including astrocytes, release purine nucleotides (the ATPs) as extracellular messengers to attract immune cells and modulate astrocytic activity⁶⁷. Imbalances in ATP levels resulting from purinergic pathway impairment have been linked to sleep deprivation and depression^68,69. As in the aforementioned pathways, ionotropic (P2 × 1-7) and metabotropic (P2Y1, 2, 4, 6, 11–14) receptors (also designated as P2 receptors), along with the other four P1 type adenosine receptors (A1AR, A2AR, A2BR, and A3AR), maintain purine nucleotide levels⁷⁰. Expression levels of P1 receptors were not altered significantly in our study (Fig. 4 panel B). Previous research has shown that HIV’s Tat protein can induce P2 × 7 (P2RX7) expression in astrocytes. This receptor is widely distributed in the CNS and is also activated during injuries and neurodegenerative diseases⁷¹. In our RNA-seq data, we found that SARS-CoV-2 infection had little downregulatory effect on P2 × 7 expression (logFC − 0.1 both variants), whereas the infection caused significant upregulation of P2 × 5 (P2RX5), particularly after infection with Omicron (logFC 0.95; Fig. 4 panel B).

A recent study showed that the spike protein of SARS-CoV-2 increased levels of metabotropic receptor transcripts P2Y1 (P2RY1), P2Y6 (P2RY6), and P2Y12 (P2RY12)⁷². However, in our investigation, the number of P2Y1 transcripts was significantly reduced following infection with the Omicron (logFC − 1.91) and Delta (logFC − 1.55) variants. P2Y6 (P2RY6) and P2Y11 (P2RY11) transcript levels were also reduced after infection of astrocytes with SARS-CoV-2, particularly after challenge with Delta (logFC − 1.01 and logFC − 0.92, respectively (Fig. 4 panel B).

Purinergic signaling is regulated by enzymes such as adenosine deaminases 1(ADA) and 2 (ADA2) which deaminate adenine in inosine, and ectonucleoside triphosphate phosphohydrolases NTPDase1 (CD39, ENTPD1) and NTPDase2 (CD39L1, ENTPD2), which hydrolyze purine and pyrimidine triphosphate and diphosphate nucleotides. Elevated ADA levels have been observed earlier in cerebrospinal fluid during neurological disorders⁷³. It has been noted recently that COVID-19 results in increased ADA2 and decreased ADA levels⁷⁴. Additionally, our data revealed that while ADA levels stayed constant, ADA2 levels were only marginally elevated (logFC 0.57 in Omicron; logFC 0.29 in Delta; Fig. 4 panel B). NTPDase1 expression also remained constant, despite the fact that astrocytic NTPDase2 transcript levels were significantly upregulated (logFC 1.08 Delta, logFC 1.06 Omicron; Fig. 4 panel B). A recent study found that after a spike protein challenge, both NTPDases were significantly upregulated⁷⁵.

Dysregulation in the neurotrophins

Expression imbalance of neurotrophins (NTs) and their receptors in astrocytes was reported during neurotropic viral infection or in virus-induced neurocognitive disorders⁷⁶. The nerve growth factor (NGF; NGF), brain-derived neurotrophic factor (BDNF; BDNF), glial cell line-derived neurotrophic factor (GDNF; GDNF), neurotrophins 3, 4, and 5 (NT-3, -4, -5; NTF3, 4, 5), neurturin (NRTN; NRTN), artemin (ARTN; ARTN), and persephin (PSPN; PSPN) are the well-known candidates in this regard. While most of these NTs showed a slight increase in the expression in infected astrocytes, NRTN expression was downregulated (logFC − 0.69 Omicron, logFC − 0.01 Delta; Fig. 5 panel A). Among NTs, the BDNF and NT-3 are being monitored frequently in CNS infection, neurodegenerative disorders and recently in patients suffering with COVID-19^{76,77,78,79,80}. In our study, expression of BDNF was elevated, although not significantly (logFC 0.4 in both variants) in infected astrocytes. In Omicron-infected astrocytes, NT-3 expression was elicited by a mere 0.79 logFC, while in Delta, its expression remained near basal levels (logFC 0.18; Fig. 5 panel A). Recent research has shown that patients with neurological impairment and long-term COVID-19 have elevated NTs levels and lower BDNF levels^81,82. In our study, NT-4 was a significantly evoked NT that was upregulated in astrocytes infected with Omicron (logFC 1.18) and downregulated in astrocytes infected with Delta (logFC − 0.67; Fig. 5 panel A). NT-4 regulates neuron survival and differentiation; its expression is ubiquitous, stable, and mostly remain unaffected by external factors. Increased levels of NT-4 in cerebrospinal fluid were recorded rarely, such as during viral meningitis associated with mumps and encephalitis associated with influenza⁸³. Thus, SARS-CoV-2-induced changes in the expression of this gene warrant further investigation.

Glial cell line-derived neurotrophic factor, the GDNF, is one of the most common neurotrophic factors produced in astrocytes in response to neurodamage⁸⁴. In our study, the expression of GDNF in infected astrocytes was found to be slightly elevated (logFC 0.71 Omicron, logFC 0.53 Delta). Similar findings were reported in a previous study, where the authors also demonstrated that GDNF expression was proportional to the severity of the infection⁸⁵. GDNF binds with the highest affinity to GDNF family receptor α-1 (GFRα; GFRA1), and it was observed that increased levels of both molecules have neuroprotective effects on astrocytes and neurons⁸⁶. Expression of GFRα1 in our study remained nearly at its basal level in infected astrocytes (logFC − 0.14 Delta, logFC − 0.07 Omicron; Fig. 5 panel A). The expression of GFRα2 (GFRA2), another receptor in the GFRα family, was markedly elevated in astrocytes challenged with the Delta variant (logFC 1.18). Interestingly, the receptor for NT-3, the TrkC (NTRK3) was also significantly upregulated in SARS-CoV-2 infected astrocytes, and its expression was evoked almost two log folds (logFC 1.93 in both variants; Fig. 5 panel A). In the recent study, the expression of TrkC was observed downregulated in neurons when infected with SARS-CoV-2 strain WA1/2020 (Wuhan variant), and upregulated when infected with Delta and Omicron⁸⁷.

Other NTs, including ARTN, PSPN, NGF, and NRTN, were also used to monitor neurological disorders, including Long COVID, in patients with COVID-19. For example, there was evidence of a significant increase in plasma persephin in post-COVID-19 patients with neurological symptoms⁸⁸, a decrease in serum NGF of post-COVID-19 patients and an increase in serum beta NGF levels during the acute phase⁸⁹, and a significant increase in ARTN plasma concentration in long-term COVID-19 patients with cognitive impairment⁹⁰. Omicron and Delta infections raised ARTN, PSPN, NGF, and NRTN expression levels in our study, but they were still below the +/- 0.9 cut-off (Fig. 5 panel A).

In addition to NTs and their receptors, astrocytes infected with the Delta variant had significantly higher levels of ARHGDIB (Rho GDP Dissociation Inhibitor Beta; logFC 0.92, logFC 0.88 Omicron) and CALML4 (Calmodulin Like 4; logFC 0.96, logFC 0.5 Omicron; Supplementary information 2 dataset 1), which are both involved in the neurotrophin signaling pathway.

Changes in the expression of proteins that form various cell-cell junctions in the blood-brain barrier and glial barrier (glia limitans)

Tight and adherens junctions of the brain microvascular endothelial cells

Astrocytes’ strategic position between neurons and endothelial cells in the neurovascular unit enables them to influence the regulation of proteins that form tight and adherens junctions (TJ/AJ), which regulates BBB formation and maintenance. Various factors released by astrocytes, such as vascular endothelial growth factor (VEGF), glial cell line-derived neurotrophic factor (GDNF, GDNF), basic fibroblast growth factor (bFGF, FGF2), and angiopoietin 1 (ANG1, ANGPT1), can shape the characteristics of the BBB^91,92. These factors promote enzymatic systems in BMECs, polarize expression of transporters, and participate in the formation of TJs⁹³. In our study, expression of VEGF (both VEGFA and VEGFB; Fig. 5 panel B) and ANG1 was slightly downregulated (logFC − 0.17 Delta, logFC 0.04 Omicron). Expression of bFGF was increased by 0.6 log fold following infection of astrocytes, whereas GDNF levels were altered as mentioned above. In addition to these, it has been demonstrated previously that a number of other molecules released by glial cells, including endothelin-1 (ET-1, EDN1), interleukins IL-1β (IL1B) and IL-6 (IL6), tumor necrosis factor (TNFα, TNF), and macrophage inflammatory protein-2 (MIP-2a, CXCL2), alter BBB permeability⁹⁴. However, the change in permeability, which results from the tight junctional opening, is temporary⁹⁴. The EDN1 (logFC 0.9 Delta, logFC 1.1 Omicron) was significantly evoked in astrocytes in our study, while IL-1β (logFC − 1.31 Delta, logFC − 1.89 Omicron) was significantly downregulated (Fig. 5 panel B). Interleukin-6 was significantly evoked by Delta (logFC 1.29), and only by 0.54 log fold when challenged with Omicron. TNFα remained unevoked in challenged astrocytes, but MIP2a, the member of chemokine with C-X-C motif ligand, was slightly upregulated by both variants (logFC 0.23 Delta, logFC 0.31 Omicron; Fig. 5 panel B).

Astrocytes are considered an indispensable element of BBB, particularly tight and adherens junctions. Earlier studies suggested that the ability of BMECs to form a BBB (mainly formed by both junctions) was not intrinsic to these cells, but that the glial cells induced this barrier property into the brain microvascular endothelium⁹⁵. This indicates the crucial role of astrocytes in the neuroinfections in maintaining the BBB integrity. Although SARS-CoV-2 is able to cross the BBB and cause CNS infection, we did not foresee the massive dysregulation of the factors in astrocytes that alter BBB permeability, such as VEGF, GDNF, bFGF and ANG1 (Fig. 5 panel B). Such massive dysregulation could be observed in other infections caused by the Zika virus or bacteria⁹⁶. It is interesting to observe in our data the overexpression of endothelin-1 (Fig. 5 panel B), which directly increases paracellular permeability and causes BBB leakage⁹⁷. On the other hand, the significant downregulation of IL-1β in astrocytes infected with Omicron or Delta is noteworthy (logFC − 1.31 Delta, logFC − 1.89 Omicron; Fig. 5 panel B), which could be beneficial to maintain BBB integrity, as IL-1β is known to downregulate TJ proteins⁹⁸. In general, it would be said that neither the Delta nor the Omicron variant used in this study affected expression of the factors that can cause drastic change in the BBB permeability at the endothelial level.

Glial barrier

As presented above, several studies have highlighted an important role of astrocytes in the control of BBB integrity; however, less is known about the second barrier beneath the endothelial cell layer, the glia limitans, which controls the CNS inflammation as well as permeability to larger molecules, especially inflammatory cells. This barrier (classically defined as membrana limitans gliae perivascularis) is formed by astrocytic endfeet and possesses astrocytic tight and adherens junctions³⁸. To this point, we were interested in expression levels of junctional proteins in infected astrocytes. The pathway enrichment analysis revealed that several members of the claudin family (CLDNs) and cingulin (CGN, CGN), related to tight junctions, had altered their expression due to SARS-CoV-2 assault (Fig. 5 panel B). Other classical tight junctional molecules, such as occludin (OCLN), claudin-5 (CLDN5), and tight junction proteins (known as Zonula Occludens, ZO-1, -2, -3, TJP1, 2, and 3) showed no significant change in their expression levels (Supplementary information 2 dataset 1). However, OCLN was downregulated in Delta-infected astrocytes (logFC − 0.56) and upregulated in the case of Omicron (logFC 0.39). Cingulin, in particular, was significantly overexpressed (logFC 1.50) when astrocytes were challenged with Omicron; however, its expression did not alter after assault with the Delta variant (logFC − 0.02; Fig. 5 panel B). Overexpression of cingulin indicates a protective response of the astrocytes against any plausible damage to astrocytic tight junctions, as cingulin and TJP1 are necessary to hold tight junction strands to cytoskeleton³⁸. It is important to note that the astrocytic tight junction is predominantly formed by claudins 1 (CLDN1) and 4 (CLDN4), which serve as structural molecules, while JAM-A (F11R) is a prominent signaling protein⁹⁹. Expression of these three candidates was not altered significantly in our study (CLDN1 logFC 0.49 Delta, logFC 0.76 Omicron; CLDN4 logFC − 0.67 Delta, logFC − 0.22 Omicron; JAM-A logFC − 0.37 Delta, logFC − 0.34 Omicron; Fig. 5 panel B). Two genes of the claudin family (CLDN16 and CLDN23) were significantly upregulated (CLDN16, logFC 1.11 Delta, logFC 0.6 Omicron; CLDN23, logFC 1.28 Delta, logFC 1.55 Omicron) in infected astrocytes (Fig. 5 panel B). Several viruses have evolved mechanisms to exploit claudins for their dissemination, as they are the key components of tight junction structure and function. Viruses either downregulate the claudins or have evolved mechanisms to promote their degradation to disseminate through several barriers (reviewed in¹⁰⁰). The significant downregulation of claudins 3 (logFC − 1.01 Delta, logFC − 0.51 Omicron), 7 (logFC − 0.3 Delta, logFC − 0.95 Omicron), and 15 (logFC − 0.44 Delta, logFC − 0.95 Omicron) in our study (Fig. 5 panel B) merits further investigation.

Pathway enrichment analysis of our data also showed that some of the genes encoding proteins associated with AJ had changed their expression significantly in a SARS-CoV-2 variant specific manner (Fig. 5 panel B), such as, CTNNA3 (catenin alpha 3, logFC − 1.85 Delta, logFC − 0.32 Omicron), PTPRB (vascular endothelial protein tyrosine phosphatase β, logFC − 0.0 Delta, logFC 1.41 Omicron), TCF7 (transcription factor 7, logFC 1.23 Delta, logFC 0.49 Omicron) and MYLPF (myosin light chain 11, logFC 0.97 Delta, logFC 0.99 Omicron), while expression of other canonical AJ components, mainly E-cadherin proteins (CDH1) that link to actin filaments via α-, β-, and γ-catenin (CTNNA1, CTNNB1, JUP, respectively), remained at basal level in infected astrocytes (Supplementary information 2 dataset 1).

Apart from the TJ and AJ proteins, we looked at the molecules that govern stabilization of the BBB. For example, the activation of Rac1 (RAC1) promotes stabilization of TJ and AJ, while RhoA (RHOA) promotes destabilization. The expression of both molecules was unchanged in infected astrocytes (Supplementary information 2 dataset 1).

Change in non-coding RNA expression

The cell uses many ways to defend itself against viral infection, including the use of non-coding RNAs (ncRNAs), which have been identified as critical regulators in antiviral immune response via the regulation of gene expressions¹⁰¹. Tens of thousands of short ( 200 nt) and long (> 200 nt) ncRNAs have been discovered so far; among them microRNAs (miRNAs), small interfering RNAs (endo-siRNAs), PIWI-interacting RNAs (piRNAs), small nuclear RNAs (snoRNAs), small tRNA-derived RNAs (tsRNAs), natural antisense transcripts (NATs), circular RNAs (circRNAs), long intergenic noncoding RNAs (lincRNAs) are prominent species^102,103. In our analysis, among the numerous identified ncRNAs, NATs (2938 of identified transcripts) and lincRNAs (3304 of identified transcripts) were the most abundant biotypes, while other identified non-coding gene biotypes were miscRNAs, snRNAs, snoRNAs including scaRNA (Small Cajal body-specific RNA), miscellaneous lncRNAs, miRNAs, etc. (supplementary information 2 dataset 8).

Changes in ncRNA expression during viral infections have been reported by a number of studies employing transcriptomic analyses^{104,105,106,107,108,109}. Changes in ncRNA expression were also investigated in various cells infected with SARS-CoV-2, as well as in patients^110,111,112, however, there is scanty data available to date describing changes in infected astrocytes. Surprisingly, our transcriptome analysis of infected astrocytes revealed that many ncRNAs had more than two-fold altered expression levels (supplementary information 2 dataset 6 and 7).

Expression changes of miRNA, lncRNA, NAT, lincRNA or circRNA may lead to neurodegenerative and neuropsychiatric disorders^113,114. This implies that virus-induced dysregulation could trigger neurocognitive disorders, as seen in COVID-19 patients. Our analysis did not detect any of the circRNA species with altered expression in infected astrocytes; however, it recorded a considerable number of dysregulated NATs and lincRNAs. Both NATs and lincRNAs are known to regulate transcriptional, post-transcriptional, and translational processes^115,116 and are associated with neurodegenerative disorders. Some of the associated NATs are BACE1-antisense (AS), BDNF-AS, GDNF-AS, HTT-AS, LOXL1-AS1, SOX21-AS or MAGI2-AS3^117,118,119. Except for HTT-AS, expressions of these NATs were altered in astrocytes in our study (Fig. 5 panel C). BACE1–AS positively regulates the expression of its counterpart, BACE1 (encodes enzyme β-secretase 1). For example, in pathological conditions increased level of BACE1-AS caused increased expression of BACE1^{120,121,122,123}. Increased expression of BACE1-AS/BACE1 was also detected in HIV-infected human astrocytes and in astrocytes infected by simian immunodeficiency virus in macaques¹²⁴. In our study, BACE1–AS/BACE1 expression was however downregulated in infected astrocytes (BACE1–AS, logFC − 0.81 Delta, logFC − 0.54 Omicron; BACE1, logFC − 0.14 Delta, logFC − 0.07 Omicron).

Expression of another NAT, BDNF-AS (antisense transcript of BDNF) was also downregulated in our experiment (logFC − 0.39 Delta, logFC − 0.44 Omicron; Fig. 5 panel C). Several studies have shown that BDNF-AS represses expression of BDNF sense transcript by changing chromatin configuration at the BDNF region, which consecutively decreases levels of endogenous BDNF protein and function^125,126,127. It was also found that downregulation or complete inhibition of BDNF-AS led to an upregulation of BDNF mRNA, subsequently increasing BDNF protein levels, suppressing neuronal cell apoptosis, and promoting neuronal outgrowth and differentiation¹²⁸. With this background, we may hypothesize that downregulation of BDNF-AS in astrocytes infected with SARS-CoV-2 may be the defense mechanism, allowing astrocytes to produce more BDNF during infection and preventing apoptosis. This hypothesis warrants further investigation.

Another abundant non-coding biotype, lincRNAs, was discovered to be upregulated in patients with elevated SARS-CoV-2 viral load. The following lincRNAs were among the most upregulated: BISPR, MIR155HG, LINC01208, U62317.2, LINC02068 and AL512306.2¹²⁹. Our findings, on the contrary, revealed that BISPR was downregulated in both variants (logFC − 0.71), MIR155HG appeared to be downregulated (logFC − 0.51) in the Delta variant but not in the Omicron variant (logFC − 0.15), and LINC01208 expression was slightly downregulated (logFC − 0.26) in the Omicron variant but not in the Delta. U62317.2 and LINC02068 expression remained unchanged (supplementary information 2 dataset 1). The only upregulated lincRNA gene was AL512306.2 (logFC 0.42 both variants) (Fig. 5 panel C). Two more LincRNA species, MALAT1 (metastasis associated lung adenocarcinoma transcript 1) and NEAT1 (nuclear paraspeckle assembly transcript 1), are recognized for their neuroprotective properties, which include enhancing cognitive function and protecting BBB function^130,131. Upregulation of MALAT1 was reported in viral infections^132,133,134, while altered expression of NEAT1 was reported in mice brain cells infected with the Hantaan virus, HIV-1, herpes simplex virus 1 or JEV^{135,136,137,138}. It was proposed recently that dysregulation of NEAT1 can cause astrocyte dysfunction and memory deficits¹³⁹. Interestingly, in our experiments both genes were evidently downregulated in the Omicron variant (MALAT1 logFC − 0.24; NEAT1 logFC − 0.38), while the Delta variant caused slight upregulation (MALAT1 logFC 0.1; NEAT1 logFC 0.06; Fig. 5 panel C).

Other ncRNAs elicited in response to infection in our study included the antisense transcript USP30-AS1, which was downregulated by both variants; NRAV, which was significantly upregulated by both variants (logFC 0.91 Delta, logFC 0.51 Omicron); and LINC-PINT RNA, which was downregulated by the Omicron variant (logFC − 0.35; Fig. 5 panel C).

3’ RNA-seq employed in this study—pros and cons

The present study and our earlier published works^{96,140,141,142,143} have effectively used 3’ RNA sequencing to comprehensively study gene expression levels in the cells challenged with pathogens. The QuantSeq 3′ mRNA-Seq (Lexogen, Austria) combined with Illumina NextSeq provides us with a minimal depth of 8 million reads per sample, which offers a cost-effective way to quantify gene expression levels, especially for large-scale studies. However, the data generated in our studies cannot capture the full complexity of alternative splicing, novel transcripts, or the tRNA fragments (tRFs) derived from tRNAs. The full RNA-seq should be used to reveal the most comprehensive view of the transcriptome; however, it is more expensive and requires more complex data analysis.