Description of study subjectsCOVID-19 patients with Asymptomatic or Symptomatic disease were recruited and their disease severity assessed based on the WHO ordinal scoring (Table 1, S Table 1)21. We included patients for whom respiratory samples had been available and therefore SARS-CoV-2 variant identification had been performed in respiratory samples collected at the time of their first COVID-19 diagnostic PCR6,22. All cases were recruited to the study within 72 h of their positive SARS-CoV-2 diagnostic test and they had no known history of COVID-19 infection. The workflow and design of the study is illustrated in Fig. 1. Transcriptional profiles were determined based on RNA arrays run on blood samples of study subjects. DEGs identified were further investigated using GO and KEGG pathways-based analyses.Table 1 Description of study population.Fig. 1Description of work flow for analysis of blood RNA transcriptional profiles in COVID-19 cases. Samples were categorized based on COVID-19 cases with Symptomatic or Asymptomatic infections. They were identified by their infecting SARS-CoV-2 variant as wildtype, or Variants of Concern (VOC). Whole blood was processed for RNA microarray analysis using the Affymetrix system and differentially regulated genes identified (DEGs) identified. These were analysed by bioinformatics analysis using KEGG GO, Reactome and WikiPathways.Of the 51 COVID-19 cases studied, the majority (80%) of cases were recruited prior to the introduction of COVID-19 vaccinations (Table 1). 31% of individuals were females. Nineteen individuals had Symptomatic whilst 32 had Asymptomatic COVID-19 at the time of PCR diagnosis. The median age of Symptomatic COVID-19 cases (53 years) was higher than those with Asymptomatic disease (39 years, p = 0.004, MWU test). 16% of those with Symptomatic COVID-19 as compared with 3% of Asymptomatic cases had history of diabetes and hypertension. C- Reactive Protein (CRP) levels were raised in those with Symptomatic disease (p = 0.004). The neutrophil lymphocyte ratio (NLR) showed a higher trend in Symptomatic patients but these data did not differ significantly likely due to the limited number of Asymptomatic cases for whom the laboratory parameters were available. SARS-CoV-2 genome sequencing of respiratory samples showed that (n = 28, 55%) patients were infected with wildtype (wt)s SARS-CoV-2. Nineteen patients were infected with VOC; four were infected with alpha, 12 with delta and three with omicron variants (other variants were 19 A and 21D). In two samples, variants could not be determined due to low viral RNA present.Differential gene expression between individuals with symptomatic versus asymptomatic COVID-19To understand host gene expression profiles associated with disease susceptibility we first compared transcriptional profiles, identifying DEGs between cases with Symptomatic (n = 19) and Asymptomatic (n = 32) COVID-19. Unhierarchical Principal Component Analysis (PCA) of transcriptional signals from the 51 COVID-19 patients showed a separation of Symptomatic COVID-19 (blue) clustered in PC1 indicating that expression profiles can differentiate health and disease in COVID-19 (Fig. 2A). Interestingly, cases with Asymptomatic infections showed equal distribution in PC2 and PC3. A total of 1443 DEGs were identified of which 608 genes were upregulated and 835 downregulated DEGs (Fig. 2B). To identify genes most differentially regulated between Symptomatic and Asymptomatic COVID-19, a volcano plot was generated (Fig. 2C) depicting genes which showed > 2-logFC change with p value < 0.05. This analysis identified the neutrophil activating chemokine CXCL8 and the IL1R2 cytokine gene the most differentially increased. In addition, there was upregulation of other gene associated with inflammatory (NFKB1A, TNFAIP3, MMP8, LTF, NFIL3, MCMEMP1), apoptotic (BCL2A1 and ANXA3) and blood coagulation (TFPI) respectively, in Symptomatic cases. Notably, there was downregulation of TCF7, a transcription factor that plays an important role in T cell development23. A heat map was generated which identified the top 30 up-DEGs (Fig. 2D) in the unclustered hierarchical map. This analysis identified differential profiles of inflammatory genes belonging to the innate cellular network involved in host defense at the initial stages of viral infection. Meanwhile, the down-DEGs related to lymphocyte activation and signaling (such as, CD3E, CD79A, CCR7, IL7R, CXCR1) as well as TCF7 involved in immune activation of T cell subsets that play an important role in controlling disease progression in the case of viral escape in the earlier phase of infection.Fig. 2Comparison of transcriptional expression between those with Symptomatic and Asymptomatic COVID-19. (A) Principal component analysis to visualize the clustering of datasets into Symptomatic, n = 19 (Symp = red), Asymptomatic, n = 32 (Asymp = blue) groups. (B) Barplots depict the 2051 differentially expressed genes (DEGs) between Symptomatic and Asymptomatic COVID-19 cases as Up-regulated, or Down-regulated and Total genes. (C) Volcano plot of DEGs between Symptomatic versus Asymptomatic COVID-19 cases. The log2 (FC) (fold change) is plotted on the x-axis, and the negative log10 p-value is plotted on the y-axis. The red points on the plot show Upregulated expression, Downregulated genes are shown in blue. Analysis with the absolute value of log2 (FC) not less than 1 (FC = 2) and p values less than 0.05. (C) The upper panel shows the 36 top-most Up-DEGs and the lower panel shows 18 Down-DEGs identified in the volcano plot. (D) Unhierarchical clustered heat-map of DEGs. (E) Dotplot of GO Analysis on biological pathways applied on Symptomatic vs. Asymptomatic COVID-19, with pathways on y-axis and gene ratio on the x-axis. (F) GSE Reactome pathway analysis applied on data. Pathways displayed on the y-axis and gene ratio on x-axis. The greater the size of circle the greater the number of genes involved in a pathway the circles are colored based on p-adjusted value. (G) A histogram of DEGs compared between Symptomatic and Asymptomatic COVID-19 cases focusing on inflammatory genes and ISGs (Table 2). The upper panel shows up-DEGs (blue) with gene IDs on the horizontal axis and fold change (log2) on the y axis. The lower panel shows down-DEGs (pink) with y axis units are fold decrease. Genes selected had < 2 > fold change with p value < 0.05 between comparative groups. (H) mRNA expression was evaluated using RT-PCR for COVID-19 cases. HuPO was used as the housekeeping gene to normalize the expression of the target genes: MAVS, OAS1, IFNa, IF144, HBB and ALAS2. Relative expression levels were quantified using the 2–∆∆CT method. The Mann-Whitney tests were performed to determine significant difference (p ≤ 0.05) between gene expression between the asymptomatic group (n = 22) and the symptomatic group (n = 4).We next conducted GO enrichment analysis of biological pathways which revealed that in those with Symptomatic COVID-19, pathways related to the antimicrobial humoral response, inflammatory response and others related cellular responses to cytokines and signal transduction were upregulated. Whilst genes related to T cell selection, activation and signaling were suppressed (Fig. 2E). This is further emphasized through the GSE reactome pathway analysis (Fig. 2F) which showed immune activation of IL-1 and IL-10 responses in addition to other cytokines identified through GSE analysis of KEGG molecular pathways.We further focused on Interferon-driven responses which play a fundamental role in determining disease outcomes in COVID-19. We investigated the presence of 96 selected genes including interferon stimulated genes (ISG), cytokine, chemokine and inflammatory pathway genes (S Table 2). By comparison of DEGs between Symptomatic and Asymptomatic COVID-19 cases we observed that in Symptomatic cases there was an upregulation of 16 genes (e.g. IL1, IL4, IL18, AREG, SOCS3, VEGFA) with downregulation of 22 genes (e.g. IL12, IL2, IFNAR2, TGF, MAVS, IRF4) compared with Asymptomatic cases (Fig. 2G).Differential gene expression of selected genes was validated through RT-PCR studies using a sub-set of RNA samples available from Symptomatic (n = 4) and Asymptomatic (n = 19) cases, Fig. 2H. mRNA levels of MAVS were higher in Asymptomatic patients (p = 0.009) whilst levels of OAS1, IFN-a and IFI44 did not significantly differ between COVID-19 groups. Expression of ALAS2 was reduced in Symptomatic patients (p = 0.02) whilst, HBB showed a lower trend but was not different. These data fit with the DEGs pattern observed in the microarray analysis (S Table 3; MAVS p = 0.0007; OAS1 p = 0.08; IFN-a p = 0.57; IFI44 p = 0.48; HBB p = 0.26; ALAS2 p = 0.059).Investigating the effect of different SARS-CoV-2 variants on the modification host transcriptional responses during COVID-19To compare how infections with SARS-CoV-2 wt compared with VOC comprising, alpha, delta and omicron strains. We studied the transcriptional profiles of 47 COVID-19 cases. First, a PCA analysis was used to view transcriptional signatures from wt (n = 28) and VOC (n = 19) infections (Fig. 3A). PCA analysis revealed clustering of individuals infected with delta and omicron, which were in a more distant compartment from profiles of those infected with wildtype and alpha variants. In total, we observed 3934 DEGs between COVID-19 cases infected with different SARS-CoV-2 variants and the analysis revealed an increasing number of DEGs with each emerging VOC; 196 DEGs between alpha and wildtype, 1425 DEGs between delta and wildtype whilst there were 2313 DEGs between omicron and wt individuals. Up DEGs were 78.6% in the case of alpha, 70.2% in the case of delta and 58.9% in the case of omicron comparisons with wt, respectively.Fig. 3Differential transcriptional expression between those infected with VOC or wildtype SARS-CoV-2. (A) Clustering of DEGs using PCA analysis of transcriptional profiles of COVID-19 cases infected with wildtype (blue, wt, n = 32), or VOC comprising; alpha (green, n = 2), delta (red, n = 11), omicron (teal, n = 2) or other (n = 4, purple) variants. (B) Venn diagram depicts overlaps between DEGs between individuals infected with wt, alpha (n = 196), delta (n = 1475) or omicron (n = 2313) strains. The numbers in the non-overlapping segments indicate unique DEGs in each comparison. (C) Comparison of ISGs genes between COVID-19 cases infected with alpha, delta and omicron variants compared with wildtype (wt), respectively. ISGs between Symptomatic (Symp) and Asymptomatic (Asymp) cases are also displayed. The red color shows the higher fold change and the lower foldchange is shown by blue color.A Venn diagram depicts the pattern of overlapping DEGs between COVID-19 VOC and wildtype infections. Only 22 genes were common DEGs between the three VOC vs. wildtype comparisons (Fig. 3B). The numbers in the non-overlapping segments indicate unique DEGs in each comparison. Despite limitations in the number of alpha and omicron, as well as non-stratification by disease category due to the small numbers in the two COVID-19 severity groups, our data suggests that the VOCs induced differential responses in the host, and we explored this further by separately comparing their gene expression against wt COVID-19 infections.Dendrogram heatmap of ISGs in the context of disease severity and infection with SARS-CoV-2 variantsActivation of type I interferon responses is associated with anti-viral immunity with increased levels associated with improved COVID-19 outcomes We compared the expression patterns of a gene signature of 47 ISGs within the DEGs identified between Symptomatic and Asymptomatic COVID-19 as well as the comparison of DEGs between alpha, delta and omicron infections, respectively. The comparative gene expression of the 47 ISGs is represented in a heatmap dendrogram (Fig. 3C). ISGs were most highly upregulated in those infected with delta, followed by omicron, alpha and then in those with Symptomatic disease. The highest fold change was observed in IF144L and IF144 in delta infected individuals, and for IRAK3 in those with Symptomatic as compared with Asymptomatic COVID-19. The lowest foldchange of the ISGs was observed for DDX5 expression in omicron infected COVID-19 cases.Alpha variants upregulate pro-inflammatory responses but downregulate protective responses as compared with wild-type SARS-CoV-2 strainsTo further investigate the 196 DEGs identified in COVID-19 patients infected with alpha (n = 4) and wt cases (n = 28), we used a volcano plot analysis. Of the 151 Up-genes, highly upregulated DEGs included SYVN1, CH25H, SAP130, VIPR1, ITM2C, TXNDC5, VPREB3, AFF3 and ITM2C, these genes may play roles in various cellular processes, including immune response, antiviral defense, and inflammation (S Table 4). SYVN1 (regulates the signal transducer and activator of transcription 3 (STAT3)/vascular endothelial growth factor (VEGF) axis) and CH25H (inflammatory factor present in macrophage populations)24 both are important inflammatory genes (Fig. 4A). Of the 42 Down-genes; CST7 (Cystatin F encoded by CST7 is a cysteine peptidase inhibitor known to be expressed in natural killer (NK) and CD8(+) T cells during steady-state conditions). GBP5 and GBP1 (guanylate binding proteins) are involved in intracellular pathogen defense mechanisms. PLSCR1 (phospholipid scramblase 1) plays roles in membrane dynamics and apoptosis. TRIM22 is an antiviral protein with roles in innate immune responses. FCGR3A (Fc-gamma receptor 3 A) is involved in immune cell activation and antibody-dependent cellular cytotoxicity. GSE analysis of biological pathways revealed that a number of innate immune and defense response related processes were suppressed in those infected with alpha variants (defense to viruses, cytokine signaling). Interestingly, no biological pathway genes were enriched in this comparison (Fig. 4B). A gene net plot of Wiki Enriched pathways further revealed that genes affected were associated with immune and inflammatory pathways with SARS-CoV-2 related signaling, apoptosis and type II interferon signaling (Fig. 4C).Fig. 4Comparison of transcriptional profiles and biological pathways in COVID-19 cases infected with alpha versus wild-type SARS-CoV-2 strains. (A) A volcano plot showing fold change in 196 DEGs between the analysis of patients infected with alpha (n = 2) versus wt (n = 32) SARS-CoV-2. The log2 (FC) (fold change) is plotted on the x-axis, and the negative log10 p-value is plotted on the y-axis. The red dot depicts up-regulated genes and blue represents downregulated genes. Highly Up- and Down- regulated genes are labelled. (B) GSE applied on GO specifically on Biological pathways on 1475 DEGs. (C) GSE dot plot analysis applied to DEGs. Or 5D. GSEA showed GO terms related to COVID-19 disease. Pathways details are displayed on the y-axis and gene ratio is plotted on x-axis. The greater the size of circle the greater the number of genes involved in a pathway; the circles are colored based on p-adjusted value. Analysis with the absolute value of log2 (FC) not less than 1 (FC = 2) and p values less than 0.05.Delta variants upregulate inflammatory and ISGs as well as viral translation as compared with wild-type SARS-CoV-2 strainsWe investigated the 1425 DEGs identified between delta (n = 12) as compared with wt (n = 28) infections. Of the 1041 Up-genes, a Volcano plot analysis revealed EGR1 (> 100 FC) to be the most Up-regulated DEG followed by IFIT (> 40 FC), S Table 5. Of the 425 Down-genes ALAS2, which catalyzes heme biosynthetic pathway (<-200 FC) to be most down-regulated DEG (Fig. 5A). A volcano plot excluding these three genes further reveals key up-DEGs as inflammatory genes (IFIT1, CCL2) and ribosomal pathway genes (RPL26, RPS24, RPL34, CCL2, RPS7, RPS21 and RPL36L and UQCRB). On the other hand, down-DEGs were CA1, MYO9B, GRN as well as erythrocyte synthesis related genes; HBB, DMTN, SCL4A1, HBD9 and GYPC. GSE analysis of GO Biological pathways revealed that the greatest number of upregulated genes were those related to cytoplasmic ribosomal proteins. Activated pathways included, VEGFA-VEGFR2 signaling, the network map of SARS-CoV-2 signaling, oxidative phosphorylation as well as nucleic acid metabolic processes (Fig. 5B). GSE dot plot analysis (Fig. 5C) which identified that all pathways were activated except for that related to multicellular organismal process, which was suppressed. Activated WikiPathways were related to RNA processing, ribosomal biogenesis, translation, peptide biosynthesis, and amide metabolic processes. These data indicate RNA, translation and ribosomal processes are activated greater in delta infected individuals, likely reflecting increased viral replication in the host. GSEA analysis of GO pathways further elaborated the enrichment score of genes involved in particular processes; with a positive regulation observed in the coronavirus disease COVID-19 pathway (Fig. 5D) as well as the GSE seen for cytoplasmic ribosomal proteins (Fig. 5E) when individuals infected with delta and SARS-CoV-2 wildtype variants were compared.Fig. 5Comparison of transcriptional profiles and biological pathways in COVID-19 cases infected with delta versus wild-type SARS-CoV-2 strains. (A) A volcano plot showing fold change in 1475 DEGs between the analysis of patients infected with delta (n = 11) versus wt (n = 32) SARS-CoV-2. The log2 (FC) (fold change) is plotted on the x-axis, and the negative log10 p-value is plotted on the y-axis. The red dot depicts up-regulated genes and blue represents downregulated genes. EGR1 was most up-regulated and ALAS2 was most down-regulated. The inset shows the volcano plot without EGR1 and IFIT1. Highly Up- and Down- regulated genes are labelled. (B) GSE applied on GO specifically on Biological pathways on DEGs. (C) GSE WikiPathway analysis applied to DEGs. GSEA showed GO terms related to COVID-19 disease. Pathways details are displayed on the y-axis and gene ratio is plotted on x-axis. The greater the size of circle the greater the number of genes involved in a pathway; the circles are colored based on p-adjusted value. Figures show the enrichment score of D, KEGG pathways that shows Positive regulation of Coronavirus disease pathways and E, WikiPathway analysis shows Positive regulation of genes involved in ribosomal protein synthesis.Omicron variants induced upregulation of viral mRNA translation and downregulation of immune responses as compared with wild-type SARS-CoV-2 strainsInvestigation of the 2313 DEGs observed between transcriptional profiles of individuals infected with omicron (n = 3) as compared with wt (n = 28) strains. A volcano plot analysis revealed EGR1 as the most Up-regulated of 1362 DEGs whilst, ALAS2 was most Down-regulated of 951 DEGs (Fig. 6A, S Table 6). Additional upregulated genes were, DNAJC15, AC4H2, MRPS21, CHD9, HMGN3 and SRP9, which suggests activation of various cellular processes, including protein folding, lipid metabolism, mitochondrial function, chromatin remodeling, and protein targeting. Of note, MGAM, a marker of innate immunity was highly down-regulated. Some of the down regulated genes were chemokine receptors (CXCR1 and CXCR2) and other inflammation related genes are TNFRSF10C, ACSL1, FPR2, SLC25A37 and ANPEP. A GSE Reactome and biological pathway analysis of the DEGs further revealed that in those infected with omicron, viral mRNA translation, cytoplasmic ribosomal pathways, gene expression, nucleic acid metabolic processes and RNA processing pathways were all activated (Fig. 6B-C). In contrast, genes involved in inflammatory and immune responses such as, chemokine signaling and viral protein interaction with cytokines were suppressed.Fig. 6Comparison of transcriptional profiles of COVID-19 cases infected with omicron versus wild-type SARS-CoV-2 strains. (A) A volcano plot of 2313 DEGs found between transcriptional profiles of cases infected with omicron (n = 2) versus wildtypes (n = 32). The x-axis is labelled with Fold change and y-axis shows –log10(p-value), the red dot depicts up-regulated genes and blue represents downregulated genes. The inset shows the volcano plot without the ALAS2 gene. (B) GSE function applied on DEGs using the Reactome pathway. (C) GSE of biological pathways. Pathways details are displayed on the y-axis and gene ratio is plotted on x-axis. The greater the size of circle the greater the number of genes involved in a pathway; the circles are colored based on p-adjusted value.Investigating VOC induced transcription in the asymptomatic hostThus far, our study of SARS-CoV-2 VOC infections showed enhanced inflammation, innate immune pathways including ISG and cytokines, with increased dysregulation of adaptive immune function. We wanted to understand how VOC may affect individuals who are able to control disease well. We therefore focused next only on studying transcriptome profiles of individuals who suffered Asymptomatic COVID-19, comprising 18 individuals infected with wt and nine cases infected with VOC (alpha, n = 4 or delta, n = 5). Of note, we excluded any vaccinated individuals in for this analysis. A PCA analysis revealed that profiles of three of the VOC cases were separated from wildtype, but two were clustered with the latter (Fig. 7A). There were 530 DEGs, with 356 Up and 174 Down DEGs between VOC and wildtype SARS-CoV-2 infections (Fig. 7B). A volcano plot analysis of the DEGs showed upregulation of AIF1, RPL26, SUB1 and SRP14. ALAS2 and SLC4A1 were highly down-regulated, as well as DMTN, GYPC, HBB, MAVS and PARP8 (Fig. 7C). GSE of biological processes identified the upregulation of cytoplasmic translation, oxidative phosphorylation and nucleotide metabolism in those infected with VOC (Fig. 7D). WikiPathway analysis of the DEGs further showed activation of ribosomal pathways, inflammatory VEGFA signaling and oxidative phosphorylation (Fig. 7E). Notably, no suppressed pathways were evident here.Fig. 7Differential transcriptional expression between Asymptomatic cases infected with VOC or wildtype SARS-CoV-2. (A) Principal component analysis helps to visualize the clustering of COVID-19 Asymptomatic cases infected with wildtype, or VOC (n = 5) strains. (B) A barplot depicts the hierarchical clustering of differentially expressed genes (DEGs) as Upregulated (red) or Downregulated (green) between individuals infected with VOC as compared with VOC. (C) A volcano plot of DEGs between infection with VOC or wt infections. The x-axis is labelled by fold change and y-axis shows –log10(p-value), the red dot depicts up-regulated genes and blue represents downregulated genes. (D) GSE of Biological processes run on the DEGs. (E) GSE WikiPathways analysis is depicted. Pathways details are displayed on the y-axis and gene ratio is plotted on x-axis. The greater the size of circle the greater the number of genes involved in a pathway; the circles are colored based on p-adjusted value.We next screened the DEGs for inflammatory pathway related genes and ISGs, observing the upregulation of inflammatory genes TNFA, CXCL8 and CCL3 but downregulation of T cell chemoattractants CXCR2 and CXCR1 in those infected with VOC as compared with wt (S Table 7). Amongst ISGs, Interferon responses were inhibited as seen by downregulation of IFIT6, IFII6, MAVS, OAS3, SAMD3, TRIM22, TRIM56. This further demonstrates that VOC infections reduced ISG responses in the host.
