Identification of PBMC-based Molecular Signature Associational with COVID-19 Disease Severity.

The longevity of COVID-19 as a global pandemic, and the devastating effects it has had on certain subsets of individuals thus far has highlighted the importance of identifying blood-based biomarkers associated with disease severity. We employed computational and transcriptome analyses of publicly available datasets from PBMCs from 126 patients with COVID-19 admitted to ICU (n=50), COVID-19 not admitted to ICU (n=50), non-COVID-19 admitted to ICU (n=16) and non-COVID-19 not admitted to ICU (n=10), and utilized the Gencode V33 assembly to analyze protein coding mRNA and long noncoding RNA (lncRNA) transcriptomes in the context of disease severity. Our data identified several aberrantly expressed mRNA and lncRNA based biomarkers associated with SARS-CoV-2 severity, which in turn significantly affected canonical, upstream, and disease functions in each group of patients. Immune, interferon, and antiviral responses were severely suppressed in COVID-19 admitted to ICU versus COVID-19 who were not admitted to ICU. Our data suggests a possible therapeutic approach for severe COVID-19 through administration of interferon therapy. Delving further into these biomarkers, roles and their implications on the onset and disease severity of COVID-19 could play a crucial role in patient stratification and identifying varied therapeutic options with diverse clinical implications.

Introduction

Since its initial reporting [1, 2] and as the COVID-19 global pandemic reaches over 150 million confirmed cases and more than 3.1 million deaths to date [3], the need for efforts into understanding the etiology of this viral infection more comprehensively with all its strains, becomes increasingly important, as well as the factors affecting disease severity. It has become apparent that severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection does not affect all individuals in the same way, ranging from asymptomatic, to a dry cough, sore throat, fever, loss of taste or smell to more severe symptoms such as difficulty breathing or shortness of breath, requiring intensive care hospital facilities, and in some cases leading to fatality resulting from severe pneumonia in 10–15% of patients, acute respiratory distress (ARDS), multi-organ failure (MOF) and intravascular coagulopathy [4].

According to early epidemiological studies, the most credible indicator of COVID-19 disease severity is age, as well as pre-existing conditions such as diabetes, hypertension and cardiovascular diseases. Such factors indicate that high risk groups are likely to possess biomarkers associated with unhealthy lifestyles and biological age, rather than chronological age [5]. To date, the molecular basis of severe COVID-19 symptoms are not well understood. Our previous studies have highlighted an important role for interferon and inflammatory response in COVID-19 [6], as well as preferential induction of chemotaxis and lipid synthesis by SARS-CoV-2 [7]. Further studies into transcriptome analysis of coding and non-coding elements, as well as convalescent blood samples of both severe and mild cases of COVID-19 should aid in identifying molecular risk factors in order to be able to predict highly susceptible individuals for severe COVID-19 infection for early intervention and tailored therapy options.

Multiple studies have identified several biomarkers which could be potentially utilized in risk stratification models for predicting severe and fatal COVID-19. Biomarkers associated with cardiac and muscle injury as well as lower enzymes were found to be significantly elevated in severe and fatal case of COVID-19 [8]. Elevated levels of IL-6 and serum ferritin in non-survivors compared to survivors in addition to elevated C-reactive protein (CRP) could partially account for the observed inflammatory response in COVID-19, that when over reactive, can lead to cytokine storms which can cause acute lung injury and sepsis, in concordance with our observation of inflammatory cytokines correlating with severe cases of COVID-19 [6, 9, 10, 11]. Other studies reported thrombocytopenia to be associated with severe COVID-19. Low levels of platelet numbers was also characteristic of higher disease severity during the SARS outbreak, occurring in up to 55% of patients with it being a significantly high risk factor for mortality [12]. Low platelet count, along with hypoxemia, have since been used by Zou et al for the development of a SARS prognostic tool with 96.2 accuracy [13]. Thirty one other observational studies involving 7613 participants, 25 of them describing the association between severe COVID-19 and low platelet counts revealing a lower platelet count in severe patients than non-severe patients. Non-survivors of COVID-19 had much lower platelet counts than those that survived [14].

Discriminatory biomarkers for diagnosis and severity prediction of COVID-19 was reported by Jiang et al showing a significant decrease in CD3+ T cells, CD4+ T cells, CD8+ T cells, and natural killer cells in severe cases admitted to the ICU compared with mild to moderate cases [15] as it has been confirmed that T cell immune responses play a vital role in recovery from SARS-Cov-2 infection, while Tregs may not play a crucial role in recovery [15, 16]. In another study by Zhao et al., RANTES (CCL5) chemokine was elevated in mild, but not severe cases of COVID-19. In contrast, IL-10 and IL-1RA were significantly upregulated in severe cases. Cytokines typically associated with cytokines storms such as IL-6 and IFN-γ were only found to be significantly elevated in later stages of severe infection, implying cytokine storms arise as the result of severe cases of COVID-19 infection rather than their cause [17]. In addition to this, another study shows 14 other cytokines to be significantly elevated in COVID-19 patients, with varying expression profiles depending on case severity. Specifically, five of these, including IP-10, MCP-3, HGF, MIG, MIP-1a, and IL-1ra were found to play a significant role in disease severity and disease progression, where levels of these markers were remarkably higher in those critically ill and with an eventual fatal outcome [18].

D-dimers concentrations can be detected in the blood post blood clot degradation, usually indicating thrombotic disorders or intravascular coagulation [19]. Several studies have since indicated D-dimer elevations in 3.75–68.0% of COVID-19 patients [10, 20, 21]. Yao et al., underwent an extensive clinical, laboratory, and radiological analysis of the characteristics of 248 cases of COVID-19 in Renmin Hospital of Wuhan University, Wuhan, China, in which they show D-dimer levels upon admission to correlate with disease severity, making it a reliable prognostic marker for in-hospital mortality [22]. Further research into the mechanisms and specificity of D-dimers in association with SARS-CoV-2 infection need to be further explored.

We previously describe the fundamental roles of long non-coding RNAs (lncRNAs) in their roles in several cellular processes. These can include the onset and progression of different diseases such as cancers and infectious disease [23, 24]. We described multiple expressed lncRNAs that were differentially expressed during the course of infection, which could serve as disease biomarkers for SARS-CoV-2 and other viral infections. Such lncRNAs including MALAT1 and NEAT1 were found to promote HIV-1 transcription and infection and promoting inflammatory mediated cell death respectively [25, 26]. Further investigation into these lncRNAs and other identified in the context of SARS-CoV-2 infection may provide valuable indication on novel methods of managing and preventing future global crises such as the pandemic we face today.

In recent weeks, many countries, particularly in Europe and the USA, have experienced a second wave of high numbers in COVID-19 infections, resulting in the spike seen worldwide. A recent study shows, and as hypothesized by the WHO, that there is no clear evidence to claim that this seasons weather or climate has a great influence on COVID-19 transmission in the most prominent urban cities in the UK [27]. On the 8th December 2020, the first Pfizer/BioNTech vaccine was administered outside of clinical trials in the UK, with millions of doses now being administered worldwide, prioritizing the elderly and health care staff [28]. With the emergence of this new vaccine, new and valuable data will be obtained on its efficiency, however this will still be limited to COVID-19, highlighting the importance of studies such as the current one for a more in depth understanding of such viral infections with the hope of preventing further global pandemics in the near future.

In our current study, we utilized computational and transcriptome analyses of publicly available PBMCs datasets from 126 patients with COVID-19 admitted to ICU (n = 50), COVID-19 not admitted to ICU (n = 50), non-COVID-19 patients admitted to ICU (n = 16) and non-COVID-19 patients not admitted to ICU (n = 10) to analyze protein coding mRNA and lncRNA transcriptome. In-depth computational analysis employing ingenuity pathway analysis (IPA) tool shed some light onto mechanistic alterations in immune portrait from severe vs mild COVID-19. We identify severely suppressed immune, interferon, and antiviral responses in COVID-19 patients admitted to ICU in comparison to COVID-19 cases which were not admitted to ICU. The identification of mRNA-based and lncRNA-based biomarkers may enable significant discrimination between milder cases and those that have significant risks of fatality. This dataset is unique in that it utilizes several subsets of patients, highlighting differences in their response, so that we may make a comprehensive comparison for novel findings. Further investigation could potentially identify several factors to be included in risk stratification models.

Materials and methods

Demographics and sample collection

The data in this study was acquired in the original study from Albany Medical Center in Albany, NY, where 126 blood samples from adult patients were acquired, as reported by Overmyer et al. [29]. Patient symptoms ranged from moderate to severe respiratory issues suspected of infection with SARS-CoV-2. One hundred patients tested positive and 26 of those tested were negative for COVID-19. Subjects were mainly of the senior population, with median age of 63. All blood samples were collected upon admission which involved the drawing quantities for two plasma preparation tubes (PPT) tubes. One tube, used for Leukocyte RNA sequencing, was processed through LeukoLOCK® filters, RNA was then eluted from LeukoLOCK filters following manufacturer recommendation, and samples were stored at -80 °C for later analyses. For more details, please refer to the original paper by Overmyera et al, [29].

Datasets and bioinformatics

Raw RNA sequencing data were retrieved from the sequence read archive (SRA) database under accession no (PRJNA660067) [29]. Detailed experimental procedures are explained in the aforementioned reference. In brief, the TruSeq Stranded mRNA kit was used for library preparation. Paired-end (2 × 50bp) sequencing was conducted on an Illumina NovaSeq6000 sequencer. Pair-end FASTQ files were retrieved using the SRA toolkit version 2.9.2 as previously described [34]. FASTQ files were subsequently mapped and aligned to the Gencode V33 assembly (both protein coding and non-coding RNAs) using KALLISTO 0.42.1 as described before [30, 31]. Normalized expression data (TPM (Transcripts Per Million) mapped reads) were sequentially imported into AltAnalyze v.2.1.3 software for differential expression analysis using 2.0-fold change and adjusted <0.05 p-value cut-off. The Benjamini-Hochberg method was used to adjust for false discovery rate (FDR), while ANOVA analysis was used for all group comparisons. Transcripts were excluded from analysis based on TPM (<1.0 raw expression value). Hierarchical clustering was performed using cosine for columns and cosine for rows, and marker finder prediction as detailed previously [30, 32].

Gene set enrichment and modeling of gene interactions networks

Ingenuity Pathways Analysis (IPA) software (Ingenuity Systems; www.ingenuity.com) was used to analyze differentially expressed genes thorough functional annotations and regulatory network analysis using upstream regulator analysis (URA). This platform analyzes upstream molecules connected to genes in the dataset through direct or indirect relationships based on changes in expression. Mechanistic networks (MN) analysis was employed by IPA to generate signaling cascades that connect upstream regulators which helps to visualize how they connect in order to explain the observed changes in gene expression. Biological processes and disease functions casually affected by the deregulation of genes in datasets can be identified via Downstream effector analysis (DEA), in addition to predictions on their activation state (Z score). IPA uses precise algorithms to predict functional regulatory networks from gene expression data and provides a significance score for each network according to the fit of the network to the set of focus genes in the database. The possibility that focus genes in the network are found together by chance is displayed as the p-value (negative log of P). The visualization of differentially affected cellular processed and pathways under various pathological conditions is made possible via the comparison analysis feature in IPA. Affected functional categories and pathways exhibiting ≥2.0 Z score (absolute) and p < 0.05 in at least one of the conditions were considered significant and were retained.

Statistical analysis

Microsoft excel 2016 and GraphPad Prism 8.0 software (GraphPad, San Diego, CA, USA) were used for statistical analyses and graphing. Two tailed t-test was used for comparative groups. Benjamini-Hochberg method was used to adjust for false discovery rate (FDR) of transcriptome data, while ANOVA analysis was used for group comparisons. Adjusted p-values ≤ 0.05 were considered significant. For IPA analyses, a Z score (−2.0 ≥ Z ≥ 2.0) was considered significant.

Results

Identification of mRNA-based biomarkers associated with SARS-CoV-2 severity

Computational and transcriptome analyses of PBMCs from 126 patients with COVID-19 admitted to ICU (n = 50), COVID-19 not admitted to ICU (n = 50), non-COVID-19 patients admitted to ICU (n = 16) and non-COVID-19 patients not admitted to ICU (n = 10) were analyzed employing the marker discovery algorithm and data retrieved from the original study by Overmyera et al., [29]. Hierarchical clustering depicted putative markers associated with each of the indicated pathological condition and enriched gene ontology (GO) associations based on mRNA expression profiles were identified. Starting with the non-COVID-19 non-ICU group, differential expression (log2) show an upregulation of GOs (upregulation; yellow) associated with immune response, defense response to viruses, regulation of I-kappaB kinase/NF-kappaB cascade, cellular response to lipopolysaccharide, antigen receptor mediated signaling pathway and transforming growth factor beta receptor signaling, while those very same categories were downregulated in all other groups (non-COVID/non-ICU, COVID-ICU, and non-COVID/ICU) (downregulation; blue. Figure 1a and supplementary table 1). Those GOs specifically upregulated in the non-COVID/non-ICU group and partially in the non-COVID-ICU group were translational elongation, viral transcription, cystolic large ribosomal subunit, translational termination and viral infectious cycle, SRP-dependent translational protein targeting to membrane and nuclear-transcribed mRNA catabolic process, and nonsense mediated decay.

Figure 1

Identification of mRNA-based biomarkers associated with SARS-CoV-2 severity. (a) Heatmap image depicting putative markers associated with each of the indicated pathological condition (COVID-19 ICU, COVID-19 no ICU, Non-COVID-19 ICU, Non-COVID-19 non-ICU) employing the marker discovery algorithm. Enriched gene ontology (GO) associations are indicated on the y axis with enrichment p values indicated. Venn diagram depicting the overlap in upregulated (b) or downregulated (c) mRNAs in CIC vs NCIC, CIC vs CNIC, CNIC vs NCNIC, and NCIC vs NCNIC). CIC: COVID-19 ICU, CNIC: COVID-19 non ICU, NCIC: Non-COVID-19 ICU, and NCNIC: Non-COVID-19 non ICU).

These categories were varied in their expression patterns in COVID/non-ICU patients and predominantly downregulated in COVID/ICU patients. GOs including ion transmembrane transport, cellular response to peptide hormone stimulus, DNA dependent negative regulation of transcription and inositol lipid-mediated signaling were suppressed in COVID/non-ICU, and even more so in non-COVID/non-ICU, whereas they were highly upregulated in the majority of COVID/ICU patients and some non-COVID/ICU patients. Some processes such as negative regulation of cell growth, and GOs related to extracellular vesicular exosome and the cytosol were upregulated in non-COVID/ICU patients and non-COVID/non-ICU patients but variably expressed, with a tendency of downregulation in the COVID-19 groups (COVID/non-ICU and COVID/ICU patients). Figure 1b and c show the overlap in upregulated (b) or downregulated (c) mRNAs in COVID/ICU versus non-COVID/ICU (CIC vs NCIC), COVID/ICU versus COVID-19/non-ICU (CIC vs CNIC), COVID-19/non-ICU versus non-COVID-19/non-ICU (CNIC vs NCNIC), and non-COVID-19/ICU versus non-COVID-19/non-ICU (NCIC vs NCNIC) in a Venn diagram. Detailed lists of our comparative mRNA findings can be viewed in supplementary tables 2–5.

Significantly affected canonical, upstream, and disease and biological function classification based on differentially expressed genes in the indicted treatment group

Based on the aberrant expression profiles of mRNAs in each of the four comparative groups (CIC vs NCIC, CIC vs CNIC, CNIC vs NCNIC, NCIC vs NCINIC, supplementary tables 2–5), ingenuity pathway analysis (IPA) highlighted multiple significantly affected canonical, upstream, and disease and biological functions where Z scores of 2.0 ≤ Z score ≤ −2.0 (depicted according to the color scale) were considered significant. Figure 2a shows multiple upregulated conical pathways (orange), particularly in the CNIC vs NCNIC group where the majority of affected pathways were activated. These include the highly upregulated cardiac hypotrophy signaling (enhanced), insulin secretion signaling pathway, integrin signaling, estrogen receptor signaling, B cell receptor signaling and many others (Figure 2a and supplementary table 6). The majority of these conical pathways are however significantly downregulated in CIC vs CNIC, including the highly downregulated pathways such as the role of Nuclear factor of activated T-cells (NFAT) in regulation of the immune response, Protein Kinase C Theta (PKCθ) signaling in T-lymphocytes, dendritic cell maturation, crosstalk between dendritic cells and natural killer cells, and calcium-induced T-lymphocyte apoptosis.

Figure 2

Significantly affected canonical, upstream, and disease and bio function classification based on differentially expressed genes in the indicted comparison groups. Enrichment heat map of canonical pathways (a), upstream regulator (b) and disease and bio function classifications (c) based on differentially expressed genes in COVID-19 ICU, COVID-19 no ICU, Non-COVID-19 ICU, Non-COVID-19 non-ICU employing the IPA algorithm. Activation Z score is depicted according to the color scale (2.0 ≤ Z score ≤ −2.0). Red indicated activation, while blue indicated suppression. Squares with a filled circle denote canonical with an absolute activation Z score <2.0.

Figure 2b (and supplementary table 7) show the predicted affected upstream regulators in each group with varying outcomes. CIC vs NCIC and CNIC vs NCNIC however appear to show similar effects on upstream regulators. As an example, both groups had activated IFNG, PRL, IFNL1, IFNA2, STAT1, TLR7 and interferon alpha, whereas both groups suppressed upstream regulators IL1RN, RC3H1, MAPK1, and IKZF1. CIC vs CNIC appears to have varied effects on upstream regulators whereas NCIC vs NVNIC appears to have no or minimal affected on upstream regulators according to IPA. CIC vs CNIC also show to suppress several disease functions including migration of cells, cell movement, leukocyte migration and immune response of cells among many other suppressed disease functions. CNIC vs NCNIC, on the other hand appears to activate the majority of disease functions such as viral infection by RNA virus, infection by lentivirus, and cell movement. Other groups show minimal effect on disease functions with an overall suppressive notion (Figure 2c and supplementary table 8).

Figure 3a, b, and c show hierarchical heat tree maps depicting affected functional categories based on differentially expressed genes, where the major boxes represent a category of diseases and functions. Figure 3 represents the COVID/ICU versus the COVID/non-ICU group and shows decreasing activation states for functions such as cellular movement, cellular growth and survival, inflammatory response, immune cell trafficking (Figure 3a), which includes suppressed activation, migration, chemotaxis, adhesion and transmigration (Figure 3b). Figure 3c shows those functions found to be activated (orange) in COVID/ICU versus the COVID/non-ICU which represent those related to infectious diseases such as replication, infectivity and viral entry.

Figure 3

Downstream effector analysis of differentially expressed genes in COVID-19 ICU vs COVID-19 non ICU. (a) Tree map (hierarchical heat map) depicting affected functional categories based on differentially expressed genes where the major boxes represent a category of diseases and functions. Each individual colored rectangle is a particular biological function or disease and the color range indicates its predicted activation state—increasing (orange) or decreasing (blue). Darker colors indicate higher absolute Z-scores. The size of the rectangles is correlated with increasing overlap significance. Illustration of suppressed immune cell trafficking (b) and activated infectious disease (c) functional categories in COVID-19 ICU vs COVID-19 non ICU.

Expression of selected markers according to COVID-19 disease severity

RNA-seq and MarkerFinder analysis also revealed several differentially expressed markers which could serve as indicators of COVID-19 severity. Those transcripts were enriched in COVID-19 patients admitted to ICU compared to all other groups include AMPH, ARG1, ARMC12, ASPH, C5orf30, CD24, CSGALNACT1, DAAM2, FKBP5, GRB10, IL1R1, IL1R2, IRAK3, IRS2, OLAH, SAP30, TPST1 and ZBTB16 (Figure 4a and supplementary table 9). Figure 4b highlights the enriched markers in patients affected by COVID-19 but have not required ICU, differentiating the severe cases from the mild ones. These markers associated with COVID-19, but have not required ICU include; CD4, CCR3, CX3CR1, FCER1A, FCRL6, ISG15, GZMA, GZMH, GNLY, PRF1, LGALS2, SIGLEC1, KLRG1, KLRK1, XAF1, IFI27, IFI44L, and USP18. All data is represented as violin lots and have ANOVA p value of <0.001 (Figure 4b).

Figure 4

Expression of selected markers according to COVID-19 disease severity. Expression of selected panel of markers enriched in the COVID-19 ICU (a) or COVID-19 non ICU (b) group. Data are presented as violin plots with the Anova p value indicated on each plot.

Identification of lncRNA-based biomarkers associated with SARS-CoV-2 severity

Our data also highlighted lncRNA based biomarkers associated with SARS-CoV-2 severity. Our findings, displayed in a heat map, depict putative lncRNA-based markers associated with each of the indicated pathological condition (COVID-19/ICU, COVID-19/non ICU, non-COVID-19/ICU, non-COVID-19/non-ICU) employing the marker discovery algorithm (Figure 5a and supplementary table 9). The heat map shows distinct clustering of markers based on each individual pathological condition with COVID/non-ICU patients showing upregulated differential expression of the first cluster of lncRNA markers (shown by pink bar on y axis), which show downregulation for all other pathological conditions. Non-COVID/non-ICU patients show upregulation for another distinct set of lncRNA markers depicted by the yellow bar on the y axis, however interestingly, these markers are also upregulated in some of the COVID cases, both those that did and did not require ICU care. Another distinct cluster of upregulated lncRNA markers was identified mainly encompassing the COVID-ICU group (blue on the y axis), where the vast majority of other pathological conditions show downregulation for this cluster of markers. The final cluster (in green Y axis) shows varied expression levels across the board, ranging from predominantly under expressed in COVID/non-ICU patients, to predominantly upregulated in non-COVID-ICU patients. Figure 5b shows the Principal component analysis (PCA) for the lncRNA transcriptome of each pathological condition. Black and green dots representing the COVID cases, both ICU and non-ICU show to mainly cluster at negative PC1 values, whereas non-COVID cases predominantly cluster towards positive PC1 values. There does seem to be however, a noticeable overlap in PC1 and PC2, indicating a degree of relatedness between the different pathological states when displayed in this manner. List of differentially expressed lncRNAs among different COVID-19 severity groups is shown in Tables 10, 11, 12, and 13.

Figure 5

Identification of lncRNA-based biomarkers associated with SARS-CoV-2 severity. (a) Heatmap image depicting putative lncRNA-based markers associated with each of the indicated pathological condition (COVID-19 ICU, COVID-19 non ICU, Non-COVID-19 ICU, Non-COVID-19 non-ICU) employing the marker discovery algorithm. (b) Principal component analysis (PCA) for the lncRNA transcriptome of each pathological condition. Expression of selected lncRNA genes enriched in COVID-19 ICU (c) or COVID-19 non ICU (d).

Upon further investigation on an individually selected lncRNA level, some lncRNA genes were significantly enriched in the COVID-19/ICU, while some were enriched in the COVID-19/non-ICU group (Figures 5c and 5d). Selected lncRNAs are presented as violin plots depicting upregulation of LINC02207, AC020636.1, AC243967.2, LINC01127, AC00727.1 and AC092746.1 in the COVID-19/ICU group, whereas those enriched in the COVID-19/non-ICU group include LINC02084, LINC02446, LINC00861, AC015911.3, LINC01871 and ANKRD44-AS1 with ANOVA p values of <0.0001.

Discussion

As the COVID-19 pandemic has surpassed an entire year of global devastation, it becomes more apparent that this will be an ongoing struggle with countless cases of progression towards severe and fatal forms of this illness. The urgent need for identification of clinical and laboratory indictors for these cases is imperative in classifying potential biomarkers for disease severity, diagnosis and therapy options. We extracted data from Overmyera et al., acquired from 126 patients with COVID-19 admitted to ICU (n = 50), COVID-19 not admitted to ICU (n = 50), non-COVID-19 patients admitted to ICU (n = 16) and non-COVID-19 patients not admitted to ICU (n = 10) [29] and went onto utilizing Gencode V33 assembly to analyze protein coding mRNA and lncRNA transcriptomes in these samples. Hierarchical clustering identified putative markers associated with each pathological condition employing the marker discovery algorithm highlighting the COVID-19 cases in need of ICU care to be enriched for gene ontologies associated with ion transmembrane transport. The general hypothesis is that ion currents play an important role in governing SARS-CoV-2 virion entry via host-membrane fusion, with Ca2+ ions being necessary in promoting fusion peptide insertion into the lipid bilayer and the endocytosis pathway [33]. Previous work on SARS-CoV and MERS-CoV show Ca2+ to be important for stimulating the fusogenic activity of the SARS-CoV fusion peptide and enhancing viral infection by approximately two-fold, respectively [34, 35]. This process entails the rearrangement of lipid bilayers further increasing Ca2+ presence [36]. In fact, studies have shown amiodarone (a Ca2+ ion channel blocker) to inhibit the spread of SARS-CoV in vitro [37]. Our recent studies have also shown preferential induction of lipid synthesis by SARS-CoV-2 [7], further promoting the notion of lipid rearrangements that could affect ion transmembrane transport, cellular responses to peptide hormone stimulus and inositol lipid-mediated signaling as indicated by our GO analysis.

For the same group (COVID-19/ICU), a significant downregulation in immune response and defense response to the virus is observed in the majority of cases, plausibly responsible for the progression into severe infection experienced by these patients. Antigen receptor mediated signaling was also suppressed, which could be owed to the apparent downregulation of the gene expression coding angiotensin-converting enzyme 2 (ACE2), the antigen facilitating the entry of SARS-CoV and SARS-CoV-2, reported in some studies, with additional cardiovascular consequences associated with ACE2 downregulation [38, 39, 40].

Comparing cases with severe symptoms requiring ICU care with COVID cases that did not require ICU, many conical pathways, upstream regulators and disease functions were affected, by either being activated or suppressed. Multiple conical pathways were found to be suppressed in ICU compared to non-ICU cases such as; natural killer (NK) cell signaling. Natural killer cells play an important role in generating an effective immune response to infection. Other studies have also highlighted NK cells to be significantly reduced in cases of severe infection compared to patients with mild infection and healthy individuals. This study suggests the use of cancer therapy drugs such as Monalizumab, as well as interferon α, chloroquine, and other antiviral agents in order to inhibit NKG2A, which restores the function of NK cells in cancers, and could also benefit those with severe cases of COVID-19 [41]. Another study also reported on depleted levels of NK cells in peripheral blood upon sudden deterioration due to COVID-19. Evidence has also shown that once patients have recovered from serious bouts of COVID-19 after therapy, NK cells numbers were restored with reduced expression of NKG2A, confirming NK cells to play a crucial role in disease severity [42]. Wilk et al., also suggests NK exhaustion to be induced by SARS-CoV-2 infection, identifying three distinct exhaustions markers on NK cells of COVID-19 patients (LAG3, PDCD1 and HAVCR2), as well as genes related to function and maturity of NK cells showing downregulated expression in COVID-19 patients [43]. Contrary to these findings, Maucourant et al. however, reported a significant increase in certain subsets of NKG2C+ NK cells in severe patients, but not in moderate COVID-19 patients [44].

While type I and type III interferon responses have previously been reported in Calu-3 host cells infected with SARS-CoV-2 in our previously published data, the most significantly downregulated upstream regulator in COVID-19/ICU patients compared to the COVID-19/non-ICU group was Interferon gamma, (IFN-γ) [6, 7]. An early study describing the clinical and immunological characteristics of 21 patients (17 male and 4 female) with COVID-19 reported the expression of IFN-γ by CD4+ T cells tended to be lower in severe cases than in moderate cases. Suppressed interferon responses in COVID-19 have shown to be a double-edged sword, with insufficient interferon responses causing immune deficiency and inability to fight the viral infection, whereas an increased interferon response and hyper-inflammation contributes to ‘cytokine storms’, leading to organ injury and sepsis [45]. Further research into the precise roles played by different interferons on their consequential responses needs to be further clarified in the context of COVID-19 disease progression and severity for the development of effective therapeutic options going forward.

In a previous study, we described the lncRNA transcriptional landscape in SARS-CoV-2 infected bronchial epithelial cells and identified a potential role for lncRNAs to be used as possible biomarkers during the course of infection. Our current study further identified aberrantly expressed lncRNAs in relation to different disease severities, where lncRNA genes were found to be enriched in the COVID-19/ICU group compared to non-ICU patients. Recently, there have been reports on the implications lncRNAs may have in the regulation of NLRP3 inflammasome and IL-6-assocciated inflammatory signaling in COVID-19, in association with several pathways like JAK/STAT, NF-κB, HIF-1α, and MAPK [46]. Nuclear enriched abundant transcript 1 (NEAT1) often reported as an oncogenic lncRNA, has been described to translocate from the nucleus to the cytoplasm to inflammasome assembly, thereby activating caspase 1 inducing inflammatory cytokine release [25]. LncRNA metastasis associated lung adenocarcinoma transcript 1 (MALAT1) has been described in several cancer types including breast [47], pancreatic [48], and lung adenocarcinoma [49], and has additionally been associated with inflammation responses and cytokine secretion, especially IL-6, promoting tissue damage [50]. Further research into the precise roles plays by different lncRNAs and how their expression affects the onset and severity of COVID-19 is needed to further clarify the precise roles played, which requires larger datasets and functional studies.

Conclusion

Expanding our understanding on what causes severe implication upon SARS-Cov-2 infection, and which factors differentiate these individuals from milder cases is highly important for more personalized, tailored prevention and treatment measures. The convenience of blood-based biomarkers associated with severity is valuable in identifying those at risk of severe symptoms leading to fatality. Our data contributes to our understanding of the internal environment in terms of differential mRNA and lncRNA expression, as well as the predicted consequences of these expression profiles with regards to affected canonical and upstream pathways related to COVID-19 infection. Further research into specific biomarkers and their implications need to be explored functionally and clinically for our continued perusal of biomarker identification and therapeutic utilization.

Declarations

Author contribution statement

Hibah Shaath: Analyzed and interpreted the data; Wrote the paper.

Nehad M. Alajez: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Funding statement

This work was supported by Qatar Biomedical Research Institute (grant number: QB13).

Data availability statement

Data associated with this study has been deposited at sequence read archive (SRA) database under accession no. PRJNA660067.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.