Bulk and Single-Cell RNA-seq Elucidates Etiology of Severe COVID-19

1. Introduction

Coronavirus Disease 2019 — COVID-19 — was a novel type of pneumonia caused by severe acute respiratory syndrome (SARS) coronavirus 2 (SARS-CoV-2). The disease rapidly emerged in Wuhan in the year 2019, and it soon spread throughout the world. This unusual disease was described one year later in 2020. In particular, its respiratory component — SARS and acute respiratory disease syndrome (ARDS) — were then characterized [1,2,3,4,5]. It was also established that the infectious agent for COVID-19 was SARS-CoV-2 — a coronavirus of probable zoonotic origin [6]. Confusingly, COVID-19 manifested an extremely wide spectrum of clinical symptoms. While most patients have mild or moderate illness, around 15% of patients progress in severity and length of the disease so much that they required hospitalization. The fatality rate in such hospitalized patients was many times higher than that of influenza [7]. Long, or severe, COVID-19 [8] was a severe and prolonged form of the disease [9]. Long COVID-19 probably led to inflammatory deterioration in function of multiple organs, including lungs, heart, gastrointestinal tract, neurological system, blood vessels, and the reproductive system. In most severe cases, patients developed ARDS, which required mechanical ventilation. Respiratory failure, sepsis, acute cardiac injury, or heart failure tended to be direct causes of death from ARDS. Thus, etiology of COVID-19 was complex, but it could be elucidated herein using unbiased genomics datasets.

Natural history of severe COVID-19 should be introduced now. The SARS-CoV-2 virus would first infect upper airways replicating in ciliated epithelial cells and inducing inflammation therein. This tended to cause sore throat, fever, and coughing. Next, the virus would strike at lungs causing pneumonia. There was severe pulmonary and systemic inflammation, which tended to cause gross damage to lung tissue. For example, alveoli would harden, inflame and fail at their role in gas exchange. Indeed, there was fibrotic scarring, lack of respiratory surfactant, and a leaky state in the lung. Moreover, coagulopathy, endotheliopathy, and vasculitis could result from the damage to the EC layer in the lungs and then spread throughout the body [10]. Thus, SARS-CoV-2 impacted negatively on many cell- and tissue-types throughout all physiological systems of the body [11]. Most importantly, the virus appeared to infect type II alveolar cells (AT2 cells) located in tiny air sacs at the ends of air tubes of the lungs. In fact, coronavirus particles could be observed in these cells using electron microscopy, and it was obvious that infected AT2 cells were damaged through senescence [12]. However, there was also ample evidence that SARS-CoV-2 attacked ECs. For example, the coronavirus could infect and damage ECs in vivo [13]. Moreover, it was shown that SARS-CoV-2 could infect and damage human umbilical vein ECs — HUVECs [14], as well as ECs in an aortic ring assay [15]. However, only single-cell genomics could provide unbiased evidence about cell types affected by SARS-CoV-2.

Based on traditional biochemical and molecular markers, it was suggested that long COVID-19 appeared related to autoimmune diseases [16]. Indeed, patients with pre-existing inflammatory or autoimmune conditions such as cancer, obesity, hypertension, cardiovascular disease, type 2 diabetes, or connective tissue disorders were at a higher risk of a fatal outcome if autoantibodies were present. One could hypothesize that the body produced high concentrations of autoantibodies to many host proteins in response to a significant viral load in many tissues. The antibodies then attacked host proteins distributed across many tissues, which resulted in a state of systemic inflammation. Indeed, presence of autoantibodies, vascular inflammation, and an inflammatory state of endothelial cells (ECs) were typical systemic symptoms of long COVID-19 [17,18]. As these changes occur over several weeks — that is times-scales suggesting transcriptional and epigenetic changes, I will review these hypotheses using unbiased genomic datasets.

Note that the genome of SARS-CoV-2 was published early in 2020 [3]. It was deposited in GenBank with a RefSeq accession number NC_045512.2. It was a single-strand ribonucleic acid (RNA) genome consisting of 29,903 base pairs. There were 11 protein-coding genes in the genome. For example, there was a multiunit replicase and a multiunit protease, and four structural proteins. Among these, protein S was a viral surface glycoprotein, which could bind to a mammalian receptor for the angiotensin-converting enzyme 2 (ACE2). Note that ACE2 could be found attached to cell membranes in the lungs, the small intestine, arteries or veins [19], but it is not generally a very abundant protein. Protein E was an envelope protein. Protein M was a membrane glycoprotein. Protein N was a phosphoprotein of the nucleocapsid. There may be additional open reading frames (ORFs), probably non-coding, which could have a regulatory role [20]. There were also at least 41 RNA modification sites on the viral RNAs [21]. Viral transcripts were typically expressed at very high levels, dominating the expression of host genes [20]. Somewhat surprisingly, the genome of the virus also encoded a furin cleavage site, which was absent from its closest known phylogenetic relatives. The furin site might enhance the tropism of the virus towards a greater number of tissue- and cell-types [22,23]. Therefore, we had sufficient genomic information about SARS-CoV-2 to investigate COVID-19 using functional genomics as early as 2020.

Herein, I aimed to review etiology of COVID-19 [24] using unbiased high-throughput functional genomics datasets. I was aided in my task by a high number and high quality of relevant genomics datasets in public databases. For example, there are hundreds of gene expression datasets relevant to COVID-19 in Gene Expression Omnibus (see Table 1). Next generation sequencing (NGS) was typically performed on messenger RNAs (mRNAs) isolated from either in vitro cell cultures or from bulk tissue samples.

2. Unbiased evidence from functional genomic studies of tissues and cells focusing on the etiology of COVID-19.

To start with, I discuss 11 whole-tissue RNA-sequencing (RNA-seq) studies focusing on the etiology of COVID-19 listed in Table 2. Bulk-sample examples were discussed below and listed in Table 2, focusing on inflammatory changes induced in peripheral blood mononuclear cells — PBMCs [25,26,27,28] during severe COVID-19, as well as within lung cells / airways [29,30,31], or in pancreatic islets [32], or within cerebrospinal fluid — CSF [33], or in conjunctival epithelium [34].

Next, I introduce 16 single-cell RNA-seq (scRNA-seq) studies, listed in Table 3, which also aim at elucidating the etiology of COVID-19 (but discriminating between individual cells). These scRNA-seq studies focused on PBMCs [35,36,37,38,39,40,41,42,43], the broncho-alveolar lavage fluid — BALF [44,45], brain [46], lungs and lung epithelium [32,47,48], heart [47], liver [47], vascular endothelium [48], or kidneys [47].

Additionally, eight meta-analyses were reviewed (see Table 4). Five of these focused on gene expression of ACE2 in scRNA-seq datasets from the lungs or airways [49,50,51,52], or kidneys [53]. Another meta-analysis built a reference library for immune cells in normal body and in different diseases [54]. Other authors looked at similarities in gene expression between COVID-19 and cancer [55]. Finally, there was also a meta-analysis of 9 scRNA-seq studies focusing on PBMCs [56].

2.1. An in vitro study in a cell line infected with SARS-CoV-2.

As early as 2020, Blanco-Melo et al. profiled transcriptomes of immortalized human cancer cell lines infected with SARS-CoV-2 and few other common respiratory viruses (such as a strain of influenza). Note that the cancer cell lines were isolated from lung tumors [29]. In short, the cells responded to SARS-CoV-2 with an inflammatory gene expression program. For example, levels of expression of chemokines were relatively high in the SARS-CoV-2 group, while levels of interferon type I and III were relatively low. Thus, these studies suggested that suboptimal transcription of antiviral countermeasures (such as interferons, inflammatory cytokines, and chemokines) played a role in the etiology of COVID-19.

In another in vitro study, Müller et al. [57] investigated islet cells from four deceased COVID-19 patients. The authors were prompted by the discovery that proteins of the virus could be detected in pancreatic tissue in the vicinity of the islets of Langerhans. Their goal was to establish whether SARS-CoV-2 could infect pancreas causing long-term damage to β-cells (potentially resulting in diabetes). An answer obtained was generally positive.

2.2. Examples of in vivo transcriptional responses in bulk tissue samples exposed to SARS-CoV-2.

The in vitro studies described above were deficient in evidence about organismal context in which a viral infection occurred. For example, status of the immune system could not be known. Neither could we know interactions between the virus and hundreds of cell types in many tissues of the body. Instead, Daamen et al. [31] followed an alternative and advantageous approach, in which they focused on in vivo activities of immune cells. Specifically, they profiled gene expression in blood samples from hospitalized patients with COVID-19. To be precise, the authors performed bulk RNA-seq on clinical samples of PBMCs. Additionally, postmortem lung tissue and postmortem samples from diseased airways were processed. Daamen et al. interpreted resulting COVID-19 datasets by discussing the importance of the expression of immunological genes. For example, many genes associated with innate immune responses were increased in expression in infected tissues. Specifically, type I interferon and other genes important for anti-viral immunity were increased in their mRNA levels. In contrast, gene expression signatures of adaptive immune response tended to be decreased.

Interestingly, not only structural cell types like epithelia, but also infiltrating immune cells were the source of transcriptional signal related to COVID-19. For example, Daamen et al. detected populations of myeloid-like cells with high inflammatory transcriptional signatures. Daamen et al. also observed that there were insufficient numbers of activated natural killer (NK) cells in diseased samples. This deficiency could prevent efficient clearing of diseased virus-laden cells. Note also that there were also insufficient numbers of regulatory, that is cluster of differentiation 8 positive (CD8+) T cells in tissues from COVID-19 patients (these cells normally mediate adaptive immunity).

2.3. Examples of RNA-seq studies of individual cells in COVID-19.

A significant advantage of single-cell analysis — scRNA-seq [58] — was in ability to identify changes in cell proportions as well as DEGs gene-by-gene and cell-by-cell. Another strategic advantage was that one could describe heterogeneity of gene expression across many cells. Moreover, viruses generally attack tissues such that infected cells are mixed with healthy ones. Recall also that there were many different cell types in each infected tissue (such as ECs, stromal cells, infiltrating immune cells, or epithelial cells). Moreover, recall that the infection was a dynamic process, and it took place in the context of co-existing pathologies (such as cancer, atherosclerosis, dementia, or chronic autoimmune disorders).

2.3.1. Expression profiling of single cells in the immune system.

Individual cells of the immune system, either circulating in the blood or present in body fluids, were analyzed to verify whether and how COVID-19 was a disease of systemic inflammation (see studies listed in Table 3). In the first of examples, Chua et al. [59] performed scRNA-seq of nasopharyngeal and bronchial samples from 19 well-characterized COVID-19 patients with either moderate or critical disease (there were also 5 healthy controls). Major epithelial cell types were identified including basal, secretory and ciliated cells, as well as FOXN4+ cells and ionocytes. Moreover, there was a subpopulation of epithelial cells with a strong interferon gamma response. There were also 13 different cell types or states of immune cells, including macrophages, dendritic cells, mast cells, neutrophils, B cells, T cells, or NK cells. In critical cases of COVID-19, there were stronger inflammatory interactions between immune and epithelial cells in comparison to moderate cases. The inflammatory interactions were identified through expression profiling of ligand–receptor pairs in both epithelial and immune cells (these were inferred from the CellPhoneDB database). These inflammatory interactions led to damage in respiratory tissues and correlated with COVID-19 severity. In other words, even more lung injury, respiratory failure, and inflammatory tissue damage was observed in critical cases. Additionally, epithelial cells from COVID-19 patients had three times greater expression of ACE2: an entry receptor for the SARS-CoV-2 virus, which correlated with interferon signalling among immune cells. In comparison to moderate cases, critical cases had inflammatory expression profiles of ligand-receptor pairs, for example inflammatory macrophages expressing many potent chemokines, chemokine ligands, interleukin 8 (IL8), interleukin -1 beta (IL1B), and tumor necrosis factor (TNF).

Similar conclusions were reached by Zhu et al. who extensively sequenced individual transcriptomes of PBMCs [60]. Five COVID-19 patients and three healthy people donated their blood for isolation of PBMCs. An atlas of single-cell gene expression was generated in both COVID-19 and influenza patients. For example, three signaling pathways were turned on in COVID-19: the apoptosis pathway, signal transducer and activator of transcription 1 (STAT1), and interferon regulatory factor 3 (IRF3) pathway. Specifically, three different mechanisms of apoptosis were more active in T cells in COVID-19: X-linked inhibitor of apoptosis (XIAP)-associated factor 1 (XAF1) pathway, TNF, and Fas receptor pathway. In influenza, signal transducer and activator of transcription 3 (STAT3) and nuclear factor kappa-light-chain enhancer of activated B cells (NF-kappa-beta) tended to be active instead of STAT1 / IRF3. According to the results of the authors, there was also an increase in the number of fully differentiated, i.e. producing specific antibodies, B cells among the PBMCs in COVID-19. There was also a decrease in the number of lymphocytes, and this effect was mediated through the apoptosis pathway. Furthermore, individual up-regulated DEGs tended to encode pro-inflammatory cytokines, cytokine receptors, or interferon-responsive transcription factors. Overall, the study suggested that the transcriptional response in immune cells during COVID-19 was both diminished in magnitude and pro-inflammatory.

Note that a similar study was performed by Wilk et al. [27] who sequenced transcriptomes of PBMCs from 7 hospitalized patients with COVID-19, four of whom had ARDS, as well as six health controls. 44,721 cells were sequenced with an average cell number of 3,194 per sample. Uniform manifold approximation and projection (UMAP) identified 30 cell clusters and DEGs were calculated for each cluster. Large phenotypic differences were identified between COVID-19 patients and controls, in cell populations of monocytes, T cells and NK cells. There were also COVID-19-related changes in cell proportions. There were also several cell types that seemed depleted in COVID-19, namely a subtype of T-cells, conventional dendritic cells (DCs), plasmacytoid dendritic cells (pDCs), CD16⁺ monocytes and NK cells. Note that in patients with ARDS only DCs, CD16⁺ monocytes and NKs were only significantly depleted. An expanding population of plasmablasts and developing neutrophils was also detected: levels of these cells were most elevated in patients with ARDS.

Another set of authors applied scRNA-seq of PBMCs to compare immune cells in the blood in even further detail. There were 57,669 high-quality transcriptomes sequenced and assembled for 7 COVID-19 patients compared against 5 healthy controls [25]. As many as 25 immune cell subsets were identified using UMAP in resulting NGS data. A few subsets of PBMCs were primarily identified in patients with COVID-19, for example monocytes and T cells with high expression of interferon-stimulated genes including interferon-alpha inducible protein 27 (IFI27), IFITM3, or ISG15. Moreover, DEGs were determined between analogous subsets in the COVID-19 group versus control. In these DEGs, there were antiviral pathways, which were induced in the COVID-19 group in monocytes and dendritic cells. Special attention was paid to the expression of interferons, however, only modest levels of interferon-gamma expression in T and NK cells were detected. Thus, a comprehensive atlas of transcriptional immune modulations in PBMCs in COVID-19 was presented. It was also confirmed that very high levels of inflammatory cytokines in COVID-19 led to severe autoimmune disease (resulting in shock, multiple organ failure, or respiratory failure).

Bernardes et al. also profiled immune responses in PBMCs, but they explored disease trajectories chronologically in a longitudinal approach [26]. The analysis included scRNA-seq data for 358,930 cells with 10,900 cells on average per sample. There were up to four longitudinal samples per patient. The goal was to identify dynamic changes in expression, which could be correlated with disease severity in COVID-19. This goal was achieved by extensive bioinformatics analyses of whole-blood scRNA-seq data along five different time points. Following sequencing, DEGs, functional trends, networks of transcription factors, and co-expression modules were identified. Generally, the observed set of DEGs suggested that there was a transcription factor activity related to inflammation and a signature of interferon signaling in COVID-19. However, interferon response seemed to be defective in critical cases. Moreover, severe COVID-19 was associated with hypoxia and hypoxic signaling (strongly increasing erythropoiesis), as well as genome-wide hypomethylation. There was an increase in the expression of genes for a pathway linked with platelet production when survivors were compared with non-survivors. There was also an induction of pro-inflammatory cytokines in the non-survivor group. Finally, increased numbers of plasmablasts, erythroid cells, and interferon-activated megakaryocytes were found to be characteristic of severe disease. When survivors were compared with non-survivors, there were 16 TFs that were differentially regulated in non-survivors only, and 7 TFs in the opposing group. Pathway analysis using the Reactome database suggested significant enrichment for the following two terms in the non-survivor group: megakaryocytes development and platelet production, and TRAF6-mediated induction of pro-inflammatory cytokines.

Another study compared cytotoxic T-cells in a severely ill group with those in a moderately affected group [61]. In the first step, RNA libraries were constructed each representing individual cytotoxic T-cells. In the second step, single-cell sequencing was performed. As a result of the sequencing, the M protein of SARS-CoV-2 was identified as a frequent target of the cytotoxic T-cell receptor (TCR). Specifically, the M198–206 amino-acid sequence was identified as an important epitope. In fact, 81.1% of the libraries, i.e. 30 out of 37, were identified as responding to this epitope. Further downstream in their analysis, Ogura et al. identified DEGs in a comparison between populations of mature cytotoxic T cells with those that were naive. For example, granzyme A and B (GZMA and GZMB) were serine proteases, which were crucial for induction of apoptosis in a cell that was attacked with cytotoxicity. Overall, the authors concluded that T cells in after severe COVID were less capable of a cytotoxic response than those in moderate COVID. However, a similar study by Meckiff et al. [39] found more cytotoxic T cells among PBMCs in hospitalized patients with COVID-19.

In another example [28], NGS sequencing was performed on single-cell transcriptomes of PBMCs from ten moderate, six critical, and five fatal cases of SARS-CoV-2. To be precise, 57,049 single-cell transcriptomes of PBMCs were sequenced. The authors identified eight functional clusters of cells among the PBMCs: (1) CD4+ T cells, (2) CD8+ T cells, (3) B cells, (4) plasma cells, (5) NKs, (6) conventional dendritic cells, (7) canonical monocytes, and (8) non-canonical monocytes. Significantly, there were greater proportions of transcriptomes derived from myeloid cells in PBMCs of fatal COVID-19 cases (in comparison with moderate cases). These myeloid cells (i.e. proliferating bone marrow cells) were characterized by biased up-regulation of a platelet-activating signature. The authors also concluded that endothelial injury and pathological thrombotic events were common in COVID-19 (and that these events were positively correlated in frequency with activity of myeloid cells). Note that two additional and related studies of PBMCs also suggested that subtle changes in proportions of myeloid cells correlated with severity of COVID-19 [62,63].

Next, Sinha et al. [64] looked at the effects of the administration of dexamethasone on circulating neutrophils during severe COVID-19. At that time, dexamethasone became a standard treatment for COVID-19, and it was given either oral or intravenous at a dose of 6 mg once daily. A total of 15 thousand single cells, which contained an equal proportion of leukocytes and lymphocytes, were analyzed. Dexamethasone expanded immunosuppressive properties of neutrophils and changed them functionally from information receivers into information providers.

In another example focusing on immune cells within the respiratory system, Bost et al. [44] performed scRNA-seq of immune cells isolated from BALF. The patients had either severe or mild COVID-19. In parallel, the authors developed a method called Viral-Track to computationally interpret their datasets. Viral-Track was a robust, unsupervised, bioinformatics method that could detect viral RNAs in large scRNA-seq datasets. The method was benchmarked on NGS datasets from virus-infected tissues. Using Viral-Track, Bost et al. compared immune cells from severe cases with those from mild cases of COVID-19. Deep impact of the virus on the immune system was detected in severe cases. For example, the authors detected marked differences in proportions of different subtypes of immune cell: either myeloid, or lymphoid, or epithelial. There was also an increase in alveolar macrophages and plasmacytoid dendritic cells (pDCs) in BALF isolated from bronchia in mild cases of COVID-19. However, tissue-resident alveolar macrophages were replaced with recruited macrophages, monocytes, or neutrophils in severe cases.

Another study by Heming et al. [33] focused on immune cells from CSF of patients who experienced heterogenous neurological symptoms during COVID-19 (this was called neuro-COVID). For example, patients with neuro-COVID could develop headache, dizziness, cognitive abnormalities, seizures, encephalitis, stroke, or brain hemorrhage. There was also an increased number of de-differentiated monocytes and exhausted cluster of differentiation 4 positive (CD4+) T cells in neuro-COVID. (Note that exhausted T cells arose because of chronic over-stimulation in a site of active inflammation.) Significantly, the exhausted T cells had diminished effector functions, and expressed co-inhibitory receptors. Finally, there was a reduced interferon response in neuro-COVID in comparison to viral encephalitis.

Another study focused on young individuals with multisystem inflammatory syndrome (MIS-C) following COVID-19 [41]. Note that MIS-C was a very dangerous inflammatory condition in which internal and external body parts became inflamed. This inflammation could be so pervasive throughout the body that it could include the lungs, heart, brain, eyes, kidneys, or the tissues of the gastrointestinal tract. Inflammatory symptoms of MIS-C could also include a cytokine storm or fever. Following scRNA-seq, transcriptomes of PBMCs were fractioned computationally using UMAP into 30 different subpopulations of immune cells. Subsequently, DEGs were identified for each of the subpopulations of PBMCs (in a comparison between cases of MIS-C and healthy controls). The goal was to identify gene markers allowing us to better diagnose and predict severity of MIS-C. For example, MIS-C patients had higher expression of S100A-family alarmins, as well as decreased gene expression signatures characteristic of antigen presentation. Finally, there was also an elevated expression level of cytotoxicity genes in NKs and CD8+ T cells.

On the other hand, Su et al. [65] investigated the etiology of post-acute sequelae of COVID-19 (PASC), using scRNA-seq and several other multi-omics technologies. This study focused on PBMCs, and it had a long follow-up period. The study found that there were four risk factors associated with the diagnosis of PASC: (1) type 2 diabetes, (2) high SARS-CoV-2 load in the plasma, (3) Epstein-Barr viremia, (4) and the presence of specific autoantibodies.

2.3.2. Expression profiling of the lung tissue.

As patients with severe COVID-19 die of pulmonary inflammation and ARDS, scRNA-seq on samples from diseased lungs were also necessary for the understanding of the etiology of COVID-19. For example, Melms et al. [30] provided a comprehensive census of lung cellular signaling and differentiation states in lethal COVID-19. They noted that individuals with COVID-19 had highly inflamed lungs with marked infiltration of activated macrophages, but also with diminished T cell responses. Impaired activities of T cells were observed, and these were likely to have contributed to fatal outcomes. Myeloid cells were also an important source of dysregulated inflammation, being more prevalent in diseased lungs in comparison to lungs from the control group. Moreover, there was significant lung fibrosis in COVID-19, and it was correlated with the length of the duration of the disease. Indeed, numerous pathological fibroblasts were observed in disorganized tissue regions. A reduction in the epithelial cell compartment was also observed, because there were fewer alveolar type I (AT1) and alveolar type II (AT2) cells in inflamed lungs, impairing their regeneration.

Another study of multiple organs led to similar conclusions about inflammation in COVID-19. Specifically, Delorey et al. generated an extensive single-cell RNA atlas from the lungs, kidneys, livers, and hearts of dead COVID-19 patients. Following transcriptional profiling, integrated computational analysis suggested that diseased lungs were characterized by extensive inflammatory damage and failed regeneration. There was evidence of such damage in cells from patients, with and without viral RNAs. Epithelial progenitors failed to regenerate sufficiently to prevent a reduction in the numbers of AT2 cells. Moreover, gene pathways that were increased in transcriptional output in COVID-19 included: apoptosis linked with oxidative stress in pericytes, various immune pathways, cell adhesion pathways, as well as the pathway for fibroblast differentiation.

2.3.3. Expression profiling of the brain, ocular epithelia, and the vasculature.

The NGS or scRNA-seq approaches were already applied in 2021 to profile brain tissues infected with SARS-CoV-2 [46]. Cells from brains of eight patients ill with COVID-19 were compared with cells obtained from 14 healthy controls. The cause of death for nearly all patients was pneumonia with extreme inflammation, which followed more than fourteen days of mechanical ventilation. The authors discovered that microglia and astrocyte subpopulations in the brain resembled those found in neurodegenerative disease. Moreover, the study was successful in linking genetic susceptibility to COVID-19 with genetic susceptibility to neurological degeneration (for example, with genes for impaired cognition, schizophrenia, or depression). In conclusion, the study found a considerable amount of evidence that long-term inflammation is a part of brain-linked COVID-19, as if inflammation from the lung tissue was relayed to the brain [46].

In another example, Jackson et al. [34] applied scRNA-seq to an in vitro model of human ocular epithelia infected with the virus. In their experiment, ocular epithelial cells were isolated from eyes donated for research from COVID-19 patients. Specifically, adult human eyes from three female donors aged 52, 78, and 80 years old were donated. The ocular epithelial cells were grown in vitro on mitotically inactivated 3T3 feeder cells. Discussing their results, Jackson et al. suggested that there was evidence for the detection of the coronavirus in the ocular epithelial cells. Moreover, conjunctival epithelium was permissive to the entry of the SARS-CoV-2 virion. However, there was no evidence of productive replication of the virus within such epithelial cells [34]. In other words, the virion of SARS-CoV-2 could probably infect more cell types than those in which it could replicate.

Another set of authors [48] characterized the vascular system in COVID-19. Before this study, ECs were already known to contribute significantly to inflammation in COVID-19, as well as to the invasion of immune cells, vascular leakage, thrombosis, and hypoxia [66]. Motivated by such observations, de Rooij et al. performed NGS on mRNAs from ECs in cases of lethal COVID-19. In the first step, single-cell transcriptomes from more than 175,000 cell nuclei were sequenced. Next, UMAP analysis of the resulting expression profiles allowed the differentiation of the nuclei into four major classes: epithelial, stromal, endothelial, or from the immune system. In the following step, a more detailed UMAP analysis identified as many as 14 different endothelial subtypes among 35,000 EC nuclei. For example, there were arterial, pulmonary vein, large vessel, systemic vein, proliferating, capillary, or lymphatic ECs. In their conclusions, the authors showed that pulmonary ECs of deceased COVID-19 patients were enriched in genes involved in cellular stress. Moreover, there were gene expression signatures in ECs that were suggestive of failed immunomodulation and impaired vessel integrity. When proportions of cells were quantified, there were more capillary and venous ECs in COVID-19 than in controls.

3. Meta-analyses of datasets related to the etiology of COVID-19.

As there are so many datasets of gene expression available for COVID-19, it was also possible and interesting to merge multiple datasets in a meta-analysis (see Table 4). Goals of a typical meta-analysis include re-analysis, recalculation of tests of significance, statistical tests of replicability and of robustness of conclusions.

In the first example, Ziegler et al. [49] focused on identifying lung cells that expressed two major receptors for the SARS-CoV-2 virus: namely ACE2 and transmembrane serine protease 2 (TMPRSS2). The authors first focused on identifying cellular targets for SARS-CoV-2 in the lungs of humans and non-human primates, then in the gastrointestinal tract, then in the upper airways. ACE2 was identified as the main target. Next, it was established that ACE2 expression was enhanced by interferons. Similarly, a meta-analysis reported by Aguiar et al. (Aguiar, Tremblay et al. 2020) focused on ACE2, and reported evidence for a mechanism dynamically regulating ACE2 expression in response to an infection with SARS-CoV-2. Another related meta-analysis looked at ACE2 levels in kidney disease associated with COVID-19 [53].

In another example, a major re-analysis and meta-analysis published by Garg et al. re-evaluated 20 published hypotheses about the nature of the immune response in COVID-19 by integrating 9 PBMC-focused datasets [27,36,56,67,68,69,70]. The cells were classified into 5 major populations: lymphoid cells, myeloid cells, B cells, epithelial cells, and platelets. However, this meta-analysis focused on the characteristics of a single cellular compartment (i.e., immune cells) in COVID-19, rather than trying to compare responses across different cell- and tissue-types. T cells tended to decrease in number with increasing severity of COVID-19. Interferon type 1 responses tended to increase in COVID-19, and B cells tended to expand in selected clones in response to antigens from SARS-CoV-2.

Furthermore, Chen et al. [55] set out to compare transcriptional responses observed in human cell lines during infection with SARS-CoV-2 with those observed in cancer samples processed by The Cancer Genome Atlas (TCGA). The following cell lines were infected with the virus: lung adenocarcinoma cell lines A549 and Calu-3, normal human bronchial epithelial cells. Many similarities were detected, suggesting that the immune cell responses and inflammatory programmes induced in host cells were similar in both diseases.

Note that smokers are strongly susceptible to severe COVID-19. To investigate the impact of smoking on COVID-19, Smith et al. [32] performed an extensive meta-analysis of expression profiles measured for lung and airway epithelial samples. Another goal of this meta-analysis was to identify cell and tissue types that have high levels of the SARS-CoV-2 receptor ACE2.

In contrast, the goal of Muus et al. [52] was to identify cell type specific associations of age, sex, or smoking, with expression levels of ACE2 and other SARS-CoV-2 receptor molecules. To this end, publicly available scRNA-seq datasets were downloaded from GEO. These datasets were searched using the following requirement: (1) that unnormalized count data were provided; (2) that data were generated using the 10X Genomics’s Chromium platform; and (3) that human samples were profiled. In effect, the authors performed a meta-analysis of 31 lung scRNA-seq studies with over 1 million 300 thousand cells, from 377 airway and lung samples, from 228 individuals. It was confirmed that cell-specific expression patterns play a role in the etiology of COVID-19.

Finally, Zhang et al. [54] set out to identify cellular phenotypes, which were shared across many inflammatory diseases including COVID-19. Through a meta-analysis and an integration pipeline, similarities in gene expression and signalling between COVID-19 and other inflammatory diseases were to be uncovered. This meta-analysis modelled also the effects of technology, tissue of origin, and donor. The meta-analysis built a reference library for > 300,000 cells in normal body and in different inflammatory diseases as well as in COVID-19. This major cross-disease study suggested that interferon-gamma and TNA-alpha signalling in macrophages were essential for the inflammatory phenotype seen in severe COVID-19. In particular, the CXCL10+ CCL2+ inflammatory macrophage state was very abundant in severe COVID-19.

4. The etiology of severe COVID-19 in the light of gene expression data.

Based on 27 datasets reviewed together with 8 meta-analyses, we arrive at mature understanding of the etiology of severe COVID-19 (see Cover Figure). Note that I considered mostly high-impact examples of RNA-seq studies of COVID-19 (the average number of citations was 372). These studies were well-interpreted through quality publications in a peer-reviewed journal with an impact factor (IF). The advantage of RNA-seq is that changes in gene expression were likely to be longer term and stable, perhaps mediated through epigenetic reprogramming, and are likely to better characterize etiology of COVID-19 than biochemical markers. Moreover, there were several inflammatory syndromes complicating long COVID-10 and described by gene expression datasets from various tissues and organs, for example neuro-COVID [33,46], ARDS, or MIS-C [41,43].

Cover Figure. Main themes of the review. In severe COVID-19, the SARS-CoV-2 virus affects the lungs and the bronchi (A). However, not only the respiratory system is affected, but also the blood, or blood vessels, or endothelial cells, the brain, or kidneys (B). Changes in gene expression in these organs can be characterized in unbiased fashion using RNA-seq or single-cell RNA-seq (C). Resulting datasets are processed using bioinformatics tools or databases and analyzed statistically (D). Based on 27 such datasets reviewed herein together with 8 meta-analyses, we arrive at mature understanding of the etiology of severe COVID-19. Created with BioRender.com.

Cover Figure. Main themes of the review. In severe COVID-19, the SARS-CoV-2 virus affects the lungs and the bronchi (A). However, not only the respiratory system is affected, but also the blood, or blood vessels, or endothelial cells, the brain, or kidneys (B). Changes in gene expression in these organs can be characterized in unbiased fashion using RNA-seq or single-cell RNA-seq (C). Resulting datasets are processed using bioinformatics tools or databases and analyzed statistically (D). Based on 27 such datasets reviewed herein together with 8 meta-analyses, we arrive at mature understanding of the etiology of severe COVID-19. Created with BioRender.com.

In general, lists of DEGs and differentially expressed pathways strongly suggested that inflammatory mechanisms were responsible for the development of severe COVID-19 (see Table 5). In particular, six pathways were differentially expressed (Table 6). There was frequently increased interferon response, interferon signaling, or interferon-responsive transcription factors (TFs). Moreover, there was a tendency for increased immune or inflammatory responses. Expression of cytokines, chemokines, or their receptors also tended to be increased. Increased interleukin-1 family signaling. There was also Increased interleukin-1 family signaling and increased hypoxic signalling. These increases were accompanied by diminished immune system regulation, diminished angiogenesis, and diminished vessel integrity.

Not only DEGs and host-viral interactions were studied, but also cellular composition of tissues and organs. For example, changes in proportions of immune or epithelial cells detected through scRNA-seq (Table 7) gave clues as to why severe COVID-19 becomes an out-of-control immune disease affecting many organs and tissue-types. Indeed, there were frequent quantitative changes in proportions of immune cells [25,31,33,44,60,61] especially among PBMCs or among infiltrating immune cells or among NKs. For example, several studies suggested that inflammatory effector cells such as antibody-producing B-cells (plasmablasts and plasma cells), or neutrophils, or cytotoxic T cells, or activated macrophages increase in proportions in severe COVID-19 (see Figure 1 and Table 7). In contrast, regulatory immune cells such as regulatory T-cells, or dendritic cells, or antigen-presenting cells decrease in proportions. These observations suggest a hypothesis according to which a powerful immune stimulation over 7-10 days at the beginning of COVID-19 gives rise to a dysregulated immune system, which attacks the lungs leading to their scaring and fibrosis (Figure 1).

Figure 1. Inflammatory changes in proportions of immune cells in COVID-19. SARS-CoV-2 replicates in the lungs in epithelial cells, macrophages, lymphocytes, or in neutrophils (A). This leads to inflammatory syndromes in various organs and tissues in case of severe COVID-19. To investigate this, PBMCs were frequently harvested from blood samples for gene expression profiling using scRNA-seq. In response to COVID-19-induced inflammation, there were changes in proportions of cell types, which resembled auto-immune diseases or immune aging (B). The following inflammatory cell types were found to increase in proportion: neutrophils, recruited monocytes, plasma cells and plasmablasts, or activated NK cells (C). In contrast, immune cells involved in modulation of immune response tended to decrease in proportion, for example regulatory T cells, or naive NK cells. Furthermore, lung’s alveoli (E) harden, inflame, and undergo scaring as a result of severe pulmonary COVID-19 (D). Created with BioRender.com.

Figure 1. Inflammatory changes in proportions of immune cells in COVID-19. SARS-CoV-2 replicates in the lungs in epithelial cells, macrophages, lymphocytes, or in neutrophils (A). This leads to inflammatory syndromes in various organs and tissues in case of severe COVID-19. To investigate this, PBMCs were frequently harvested from blood samples for gene expression profiling using scRNA-seq. In response to COVID-19-induced inflammation, there were changes in proportions of cell types, which resembled auto-immune diseases or immune aging (B). The following inflammatory cell types were found to increase in proportion: neutrophils, recruited monocytes, plasma cells and plasmablasts, or activated NK cells (C). In contrast, immune cells involved in modulation of immune response tended to decrease in proportion, for example regulatory T cells, or naive NK cells. Furthermore, lung’s alveoli (E) harden, inflame, and undergo scaring as a result of severe pulmonary COVID-19 (D). Created with BioRender.com.

Finally, there was also transcriptional evidence for endothelial injury and pathological thrombotic events [28,48].

Note also that the transcription of the SARS-CoV-2 virus itself could be verified using RNA-seq and specialized bioinformatics pipelines, for example in epithelial cells and in lymphocytes, or macrophages, or neutrophils in bronchi or lungs of COVID-19 samples (see Table 8).

5. Conclusions.

A mean time from diagnosis to death in severe COVID-19 was estimated to be between 17 and 19 days [71], suggesting sufficient time for changes in gene expression or epigenetic regulation. Thus, severe COVID-19 is not only a viral infection with SARS-CoV-2, but also a weeks-long inflammatory transition occurring throughout the body. This transition is accompanied by global changes in gene expression patterns, particularly among PBMCs. Many tissues also undergo inflammatory transformation in cellular composition, which is accompanied by fibrosis or loss of cells that can regulate immunological reactions.