Distinct immune responses in patients infected with influenza or SARS-CoV-2, and in COVID-19 survivors, characterised by transcriptomic and cellular abundance differences in blood

Ethics and consent

The study was approved by the South Central - Hampshire A Research Ethics Committee (REC): REC reference 20/SC/0138 (March 16^th, 2020) for the COVID-19 point of care (CoV-19POC) trial; and REC reference 17/SC/0368 (September 7^th, 2017) for the FluPOC trial. For full inclusion and exclusion criteria details see (25) and (26). Patients gave written informed consent or consultee assent was obtained where patients were unable to give consent. Demographic and clinical data were collected at enrolment and outcome data from case note and electronic systems. ALEA and BC data management platforms were used for data capture and management.

Study design and participants

All participants were recruited within the first 24 hours of admission to two large studies of molecular point-of-care testing (mPOCT) for respiratory viruses (CoV-19POC and FluPOC). Blood samples including whole blood in PAXgene Blood RNA tubes (BRT) (Preanalytix) were collected from SARS-CoV-2 positive patients and influenza positive patients, within 24 hours of enrolment, and stored at −80°C. The studies were prospectively registered with the ISRCTN trial registry: ISRCTN14966673 (COV-19POC) (March 18^th, 2020), and ISRCTN17197293 (FluPOC) (November 13^th, 2017). The COV-19POC study was a non-randomised interventional trial evaluating the clinical impact of mPOCT for SARS-CoV-2 in adult patients presenting to hospital with suspected COVID-19, using the QIAGEN QIAstat-Dx PCR testing platform with the QIAstat-Dx Respiratory SARS-CoV-2 Panel (27). The trial took place during the first wave of the pandemic, from 20th March to 29th April 2020, and prior to the availability of approved COVID-19 treatments. All patients were recruited from the Acute Medical Unit (AMU), Emergency Department (ED) or other acute areas of Southampton General Hospital. The FluPOC study was a multicentre randomised controlled trial evaluating the clinical impact of mPOCT for influenza in hospitalised adult patients with acute respiratory illness, during influenza season, using the BioFire FilmArray platform with the Respiratory Panel 2.1 (28). The trial took place during influenza seasons over the two winters of 2017/18 and 2018/19. All patients were recruited from the AMU and ED of Southampton General Hospital and Royal Hampshire County Hospital.

Extraction of RNA from clinical samples and Illumina sequencing

Total RNA was extracted from PAXgene BRT using the PAXgene Blood RNA Kit (PreAnalytix), according to the manufacturer’s protocol at Containment Level 3 in a Tripass Class I hood. Extracted RNA was stored at −80°C until further use. Following the manufacturer’s protocols, total RNA was used as input material into the QIAseq FastSelect–rRNA/Globin Kit (Qiagen) protocol to remove cytoplasmic and mitochondrial rRNA and globin mRNA with a fragmentation time of 7 or 15 minutes. Subsequently the NEBNext® Ultra™ II Directional RNA Library Prep Kit for Illumina® (New England Biolabs) was used to generate the RNA libraries, followed by 11 or 13 cycles of amplification and purification using AMPure XP beads. Each library was quantified using Qubit and the size distribution assessed using the Agilent 2100 Bioanalyser and the final libraries were pooled in equimolar ratios. Libraries were sequenced using 150 bp paired-end reads on by an Illumina® NovaSeq 6000 (Illumina®, San Diego, USA). Raw fastq files were trimmed to remove Illumina adapter sequences using Cutadapt v1.2.1 (29). The option “−O 3” was set, so that the 3’ end of any reads which matched the adapter sequence with greater than 3 bp was trimmed off. The reads were further trimmed to remove low quality bases, using Sickle v1.200 (30) with a minimum window quality score of 20. After trimming, reads shorter than 10 bp were removed.

Data QC and alignment

Quality control (QC) of read data was performed using FastQC (31) (v0.11.9) and compiled and visualised with MultiQC (32) (v1.5). Samples with <20 million total reads were excluded from further analysis. The STAR index was created with STAR’s (33) (v2.7.6a) genomeGenerate function using GRCh38.primary_assembly.genome.fa and gencode.v34.annotation.gtf (34) (both downloaded from GENCODE), with –sjdbOverhang 149 and all other settings as default. Individual fastq files were aligned using the --twopassMode Basic flag, with the following parameters specified (following ENCODE standard options): --outSAMmapqUnique 60, outFilterType BySJout, -- outFilterMultimapNmax 20, --alignSJoverhangMin 8, --outFilterMismatchNmax 999, -- outFilterMismatchNoverReadLmax 0.04, --alignIntronMin 20, --alignIntronMax 1000000, -- alignMatesGapMax 1000000 and all other options as default. For rMATs (35) (v4.1.0) analysis, STAR was run again as before, but with the addition of --alignEndsType EndToEnd. Samtools (36) (v1.8) was used to sort and index the aligned data.

Comparisons of baseline clinical characteristics

Baseline clinical characteristics of the patient groups were assessed using R (37) (v4.0.2) and RStudio (38) (v1.3.959) for comparisons between COVID-19 versus influenza, and COVID-19 survivors versus non-survivors. Extreme outliers (values < Q1 - 3 interquartile range, or > Q3 + 3 interquartile range) were identified with the R package rstatix (39) (v0.7.0) and removed. Statistical testing was performed including a Shapiro-Wilk test to assess for data normality followed with either an unpaired parametric T-test (Shapiro-Wilk test p-value > 0.05) or an unpaired non-parametric Wilcoxon test (Shapiro-Wilk test p-value < 0.05) for continuous data, or a Chi-square test for categorical data. The R package Table1 (40) (v1.3) was used to plot the baseline clinical characteristics.

Systems immunology-based analysis of blood transcript modules

Blood transcript module (BTM) analysis was performed with molecular signatures derived from 5 vaccine trials (41) as a reference dataset, and BTM activity was calculated using the BTM package (41) (v1.015) in Python (42) (v3.7.2) using the normalized counts as input. Module enrichment significance was calculated using CAMERA (43) (v3.46.0). The significance threshold for the linear model was set at false discovery rate (FDR) 0.05 for the comparison between patients with COVID-19 or influenza.

Unbiased gene clustering analysis

Gene co-expression analysis was performed with BioLayout (44) (v3.4) using a correlation value of 0.95, other settings were kept at default. Clusters were manually assessed to determine gene expression differences depending on for example patient cohort. Genes were subsequently analysed with ToppGene (45) gene list enrichment analysis using the default settings.

Differential gene expression analysis between patient groups

HTSeq (46) (v0.11.2) count was used to assign counts to RNA-seq reads in the Samtools sorted BAM file using GENCODE v34 annotation. Parameters used for HTSeq were --format=bam, --order=pos, -- stranded=reverse, --type=exon and the other options were kept at default. EdgeR (47) (v3.30.3) was used for differential gene expression analysis with R (v4.0.2) in RStudio (v1.3.959). Genes with low counts across all libraries were filtered out using the filterByExpr command. Filtered gene counts were normalised using the trimmed mean of M-values (TMM) method. Differentially expressed genes were identified, after fitting the negative binomial models and obtaining dispersion estimates, using the exact test and using a threshold criteria of FDR p-value < 0.05 and log2 fold change < −1 and > 1. Genes which were within the threshold criteria were used for ToppGene gene list enrichment analysis. A principal component analysis (PCA) graph was constructed based on all differentially expressed genes to assess sample clustering.

Assessment of difference in adaptive immune response related gene expression

A higher abundance of transcripts from 83 immunoglobulin genes, overlapping with the genes in the Gene Ontology (GO) (48, 49) biological process term ‘adaptive immune response’ (Additional file 1), was found in patients with COVID-19 compared to influenza. To assess gene transcript abundance differences for these 83 genes in each patient a heatmap was generated and Z-scores were summed to give an overall positive (high) or negative (low) total Z-score. Patient baseline clinical characteristics were explored, as above, for any explanatory factors for the involvement of a high or low total Z-score between patients with COVID-19 or influenza, and those that survived COVID-19 versus those that died within 30 days of hospital admission. Metadata comparison plots were made with the R package ggplot2 (50) (3.3.2) and statistical testing with the R package ggpubr (51) (v0.4.0).

Topological mapping of global gene patterns

TopMD Pathway Analysis (52) was conducted using the differential transcript abundances identified by differential gene expression analysis, generating a map of the differentially activated pathways between all patients with COVID-19 or influenza. The TopMD pathway algorithm measures the geometrical and topological properties of global differential gene expression embedded on a gene interaction network (53). This enables plotting and measurement of the differentially activated pathways through extrapolation of groups of mechanistically related genes, called TopMD pathways. TopMD pathways possess a natural hierarchical structure and can be analysed for enriched GO terms, by chi-square test.

Assessment of differential splicing between patient groups

Three different tools were used to assess differential gene splicing between patients with COVID-19 or influenza, and COVID-19 survivors or non-survivors after 30 days of hospital admission. rMATs (35) (v4.1.0) was run using BAM files with soft clipping suppressed, generated with STAR and GENCODE v34 gene annotation. Additional settings used were -t paired, --readLength 150 and --libType fr-firststrand. Results were filtered for FDR p-value < 0.05. LeafCutter (54) (v0.2.9) was run in stages following the Differential Splicing protocol (55) (bam2junc.sh generated junction files from BAMs, leafcutter_cluster.py grouped junctions into clusters, leafcutter_ds.R tested for differential splicing, all with default settings, except –min_samples_per_intron was set to be approximately 60% of the smaller group size for each comparison (46 for COVID-19 vs influenza, 9 for COVID-19 survivors vs non-survivors), and results were filtered to exclude events with delta PSI <10%, based on recommendations (56). The LeafViz script (57), prepare_results.R was used to generate a data table from which gene names for significant events were extracted, while the map_clusters_to_genes R function was used to assign genes to non-significant tested events. Overlap between LeafCutter differentially spliced and EdgeR differentially expressed genes was tested for significance using Fisher’s Exact Test (fisher.test in R (v3.5.1) using a 2×2 contingency table and two.sided alternative hypothesis). MAJIQ (58) (v2.2) was run in two stages (majiq build and majiq deltapsi) with default settings, and results were filtered (delta PSI >20%, probability >0.95) using Voila (58) (v2.0).

In silico immune profiling predicting immune cell levels between patient groups

Relative abundance of 22 immune cell types and their statistical significance was deconvoluted from whole blood using the reference gene signature matrix (LM22) using CIBERSORTx (59). CIBERSORTx analysis was conducted on the CIBERSORTx website (60) using 100 permutations. Immune cell distribution between the groups were compared by Mann–Whitney test.

Identification of immune signatures as a predictor for COVID-19 outcome

Transcript to transcript gene co-expression network analysis with BioLayout 3D (v3.4) (Pearson coefficient 0.85, MCL=1.7) assembled 537 genes differentially expressed (EdgeR, FDR < 0.5 and |log2 fold change > 1|) in blood taken on admission between patients with COVID-19 who either survived or died of COVID-19 within 30 days of admission to hospital. Combinations of 100 genes from the top 4 clusters were assessed as predictor variables for outcome using Boosted Logistic Regression, Bayesian Generalised Linear and RandomForest models within SIMON (61) (v0.2.1) installed with Docker (62) (v20.10.2). TMM normalised gene expression data was centred and scaled. Covariant features were removed based on correlation analysis. Samples were randomly split into train:test subsets at the ratio 75%:25%.

Number of participants

In total RNA-seq was done for 80 patients with COVID-19 and 88 patients with influenza. Five patients with influenza failed QC (read count < 20M) leaving 83 patients with influenza for analysis, of which 76% were infected with the influenza A virus and 22% with influenza B virus. Two patients with COVID-19 were identified by PCA as outliers, subsequent assessment revealed an elevated white blood cell and lymphocyte count caused by pre-existing underlying chronic lymphocytic leukaemia, and these patients were excluded from further analyses (Supplementary figure 1). This left 78 patients with COVID-19, of whom 62 survived and 16 died within 30 days of hospital admission, and 83 patients with influenza.

Clinical differences

The baseline clinical characteristics of the patients used in this study for the comparison between influenza and COVID-19 were assessed and no differences in distribution of sex or age were detected between patient groups, however, more patients with influenza were of White British ethnicity (p-value 1.12×10⁻⁰⁵) and more were current smokers (p-value 9.07×10⁻⁰⁵). There were also differences in the proportion of cases with underlying comorbidities, with patients with COVID-19 more commonly having hypertension (p-value 1.42×10⁻⁰²), liver disease (p-value 3.63×10⁻⁰²) and diabetes mellitus (p-value 6.44×10⁻⁰³) than those with influenza. However, underlying respiratory disease was more common in patients with influenza (p-value 1.22×10⁻⁰³). Patients with COVID-19 generally exhibited more severe clinical symptoms. While the National Early Warning Score 2 (NEWS2) was not different between patients with COVID-19 or influenza, patients with COVID-19 had a higher respiratory rate (p-value 2.79×10⁻⁰²) and a greater proportion of patients with COVID-19 were on supplementary oxygen at hospital admission (p-value 6.81×10⁻⁰³). Laboratory results indicated higher levels of C-reactive protein (p-value 1.73×10⁻⁰³) and lymphocytes (p-value 2.76×10⁻⁰²) in patients with COVID-19. Furthermore, COVID-19 patients had a longer duration of symptoms prior to presentation to hospital (p-value 1.17×10⁻⁰⁵) and once admitted a longer length of stay (p-value 5.51×10⁻¹⁰). Longer stay time was associated with increased 30 day mortality after hospital admission and patients with COVID-19 were more likely to have died compared to patients with influenza (p-value 4.42×10⁻⁰⁵) (Table 1).

Between patients with COVID-19 who survived and those who died, a fatal outcome occurred in older patients (p-value 2.58×10⁻⁰⁹). COVID-19 non-survivors also had a shorter duration of symptoms before being admitted to hospital (p-value 5.38×10⁻⁰³). COVID-19 non-survivors more commonly had underlying comorbidities including hypertension (p-value 1.93×10⁻⁰³), cardiovascular disease (p-value 3.97×10⁰³), diabetes mellitus (p-value 2.31−10⁻⁰²) and underlying respiratory disease (p-value 1.06×10⁻⁰²). While the NEWS2 scores were not different, COVID-19 survivors had a higher heart rates than COVID-19 non-survivors (p-value 9.27×10⁻⁰³). Laboratory results showed an increase of white blood cell count (p-value 3.83×10⁻⁰²), total protein levels (p-value 2.5×10⁻⁰³), creatinine (p-value 3.87×10⁻⁰²), alanine aminotransferase levels (p-value 2.85×10⁻⁰²), troponin levels (p-value 2.37×10⁻⁰⁴), tumour necrosis factor α (TNFα) (p-value 1.43×10⁻⁰²), interleukin (IL)-6 levels (p-value 2.78−10⁻⁰³), IL-8 (p-value 2.24×10⁻⁰²), IL-1β (p-value 3.78×10⁻⁰²) and IL-10 (p-value 7.51×10⁻⁰²) in COVID-19 non-survivors. Patient outcome and length of hospital stay were different due to separation based on patient survival (Table 2).

Molecular differences

RNA-seq was used to investigate potential blood transcriptomic signatures of immune activation between patients infected with SARS-CoV-2 versus influenza, and COVID-19 survivors versus those that died. A median of 60.4 million reads in patients with COVID-19, and 58.9 million reads in patients with influenza was obtained (Supplementary figure 2A). In patients who died of COVID-19 a median of 55.7 million reads was obtained and for COVID-19 survivors the median was 62.6 million reads (Supplementary figure 2B). Clustering analysis between patients with COVID-19 or influenza indicated a homogeneity of blood transcriptome profiles suggesting any variation between groups to be subtle (Supplementary figure 3A). A partial separation was found between patients who survived or died of COVID-19 based on patient outcome after 30 days of hospital admission, indicative of a larger variation in the blood transcriptome (Supplementary figure 3B).

Contrasting innate and adaptive immune programmes

Previous studies have suggested that severe COVID-19 is associated with aberrant immune pathology (63, 64), and therefore BTM analysis and gene co-expression analysis were used to investigate the balance between the innate and adaptive response in patients with either COVID-19 or influenza virus and to identify patterns of changes associated with each arm of the immune system. A systems immunology-based analysis of BTMs between patients with COVID-19 or influenza revealed several differences (Figure 1). For the upregulated BTMs in COVID-19, signatures were observed related to the cell cycle and adaptive immune response, primarily CD4+ T cells, B cells, plasma cells and immunoglobulins. In contrast, the downregulated BTMs showed signatures associated with monocytes, inflammatory signalling and an innate antiviral and type I IFN response.

• Download figure
• Open in new tab

Figure 1: Blood transcript module (BTM) analysis between patients with COVID-19 or influenza.

Upregulated signatures in COVID-19 patients are associated with cell cycle and an adaptive immune response, primarily CD4+ T cells, B cells, plasma cells and immunoglobulins. While the downregulated signatures, associated with influenza patients, are involved with monocytes, inflammatory signalling and an innate antiviral and type I interferon response.

Gene co-expression analysis was done on a total of 4,093 transcript abundances for unbiased gene clustering between patients with COVID-19 or influenza (TMM normalised counts per million EdgeR FDR p-value <0.05) and identified 50 clusters of 4 or more genes (BioLayout 3D, Pearson R >= 0.85, MCL = 1.7). These clusters are clearly separated into groups comprising increased transcript abundances in blood of patients with influenza or COVID-19 (Figure 2) and the top 12 clusters are shown in Table 3.

• Download figure
• Open in new tab

Figure 2: Top 12 clusters identified with BioLayout.

A) Enrichment of gene clusters in blood of patients with influenza (annotated in red) and COVID-19 (annotated in blue). Increased abundances of gene transcripts in influenza patients are involved with an innate immune response, while in COVID-19 clusters are involved with an adaptive immune response, blood coagulation and neutrophil degranulation. B) After TMM normalisation a significant difference in gene clusters between patients with influenza or COVID-19 was detected. The abundance of gene transcripts involved with an innate immune response and plasmacytoid dendritic cell were observed to be higher in influenza patients. In contrast, the abundance of gene transcripts involved with an adaptive immune response and neutrophil degranulation was higher in COVID-19 patients.

Interestingly, an increased abundance of gene transcripts in patients with COVID-19 are involved in adaptive immunity, pointing to activation/priming of T cells and B cells, including induction of proliferation (cluster 4, FDR p-value 3.97×10⁻⁵⁷), Additionally, an increased abundance of gene transcripts encoding neutrophil degranulation (cluster 9) and blood coagulation (cluster 6) clearly differentiated patients with COVID-19 from patients with influenza (FDR p-value 4.33×10⁻¹⁹ and FDR p-value 2.84×10⁻¹² respectively). In contrast, an decreased abundance of gene transcripts in the blood transcriptome of patients with COVID-19 in comparison to patients with influenza were associated with innate immunity, including biological processes involved with defence response to virus (cluster 2) (FDR p-value 1.34×10⁻³⁷), type 1 helper T cell stimulation (cluster 10) (FDR p-value 4.53×10⁻⁰³), dendritic cell morphogenesis (cluster 11) (FDR p-value 1.37×10⁻⁰²), and myeloid cell activation (clusters 1 and 8) (FDR p-value 5.16×10⁻¹³ and FDR p-value 4.15×10⁻⁰⁴ respectively). Importantly, the largest decrease of transcript abundances in patients with COVID-19 comprised genes expressed in plasmacytoid dendritic cells (pDC) (FDR p-value 4.17×10⁻²²), indicating impaired immune responses to viruses (FDR p-value 1.34×10⁻³⁷) and impaired IFN signalling (FDR p-value 5.56×10⁻³⁰). This was suggestive of contrasting innate and adaptive immune programmes between the different infections and these were further investigated.

High abundance of immunoglobulin genes associated COVID-19

A total of 20,542 abundance measures of gene transcripts were obtained after filtering out transcripts with low counts, of which 4,094 transcripts were found to be significantly different between patients with COVID-19 or influenza (FDR p-value < 0.05) of which, 197 transcripts exceeded a log2 fold change of < −1 or >1, with 126 transcripts showing higher abundance in patients with COVID-19 and 71 transcripts showing higher abundance in patients with influenza (Figure 3A and Additional file 2). Complimentary to the findings from gene co-expression analysis, the transcripts with increased abundance in patients with COVID-19 were found to be involved with humoral immune response, complement activation and B cell mediated immunity (Figure 3B).

• Download figure
• Open in new tab

Figure 3: Increased adaptive immune response in patients with COVID-19 compared to influenza.

A) Volcano plot of transcripts falling within the threshold values (FDR < 0.05 and log2 fold change < −1 or >1) which were used for enrichment analysis with ToppGene. B) Enrichment analysis of the transcripts with an increased abundance in patients with COVID-19 identified an increased adaptive immune response. Percentage in annotation is the ratio of the input query genes overlapping with the genes in the pathway database. C) Heatmap of 83 immunoglobulin gene transcripts, which overlap with the GO biological process term ‘adaptive immune response’, found at a higher abundance in patients with COVID-19. Positive Z-scores are seen mostly in patients with COVID-19 while negative Z-scores are mostly seen in patients with influenza.

83 immunoglobulin genes, associated with the GO biological process term ‘adaptive immune response’, were found to have higher transcript abundance in the majority of patients with COVID-19 than those with influenza (p-value < 2.22×10⁻¹⁶, Wilcoxon test) (Figure 3C and Supplementary figure 4) and by using a total Z-score, patients with COVID-19 or influenza were classified as having either a high or low abundance of these 83 immunoglobulin genes. A high abundance was associated with a total positive Z-score (1.46 to 175.46) which was identified in 59 patients with COVID-19 and 21 patients with influenza indicating a higher than average abundance of these 83 adaptive immune response related immunoglobulin genes. While a low abundance was associated with a total negative Z-score (−0.12 to −154.93) identified in 19 patients with COVID-19 and 62 patients with influenza indicating a lower than average abundance of adaptive immune response related immunoglobulin genes. COVID-19 patients with lower abundance of adaptive immune response related immunoglobulin genes, a total negative Z-score, were found to be significantly older (p-value 6.32×10⁻³, T-test) and had a shorter duration of symptoms before being admitted into hospital (p-value 5.9×10⁻⁰⁴, Wilcoxon test). Additionally, COVID-19 patients with high abundance of adaptive immune response related immunoglobulin genes, a total positive Z-score, were significantly more likely to be still alive 30 days after admitted into hospital (x² 13.39 and p-value 2.52×10⁻⁰⁴, Chi-square test) (Figure 4).

• Download figure
• Open in new tab

Figure 4: Comparison of 83 immunoglobulin gene abundances in patients with COVID-19 or influenza.

Metadata of patients with COVID-19 or influenza with a low or high abundance of the 83 immunoglobulin genes related to GO biological process term ‘adaptive immune response’ were compared. Significant differences were detected based on age (patients with COVID-19, p-value 1.8×10⁻³, Wilcoxon-test) and duration of symptoms (patients with COVID-19, p-value 5.9×10⁻⁰⁴ and patients with influenza, 1.3×10⁻⁰³, Wilcoxon test), with older individuals and shorter symptom duration associated with the low immunoglobulin gene abundance group for COVID-19. Additionally, a low abundance of immunoglobulin genes was associated with decreased COVID-19 survival (x² 13.39 and p-value 2.52×10⁻⁰⁴, Chi-square test).

Topological mapping of global gene patterns

Topological analysis allows the measurement of the global profiles of transcript abundances relative to gene pathways without data reduction and this was used to define a global map of differentially activated pathways between COVID-19 and influenza. The first differentially activated TopMD pathway was enriched for ribosomal and insulin related pathways, with peak gene UBA52: named by GO analysis as cytoplasmic ribosomal proteins (adjusted p-value 1.55×10⁻¹⁴⁶). This pathway was also found to be enriched for genes expressed by transcription factor Myc (adjusted p-value 7.07×10⁻⁵³) against the ChEA 2016 transcription factor database and of dendritic cells in the ARCHS4 transcription factors’ co-expression database (adjusted p-value 1.34×10⁻³⁶). Activated Myc represses interferon regulatory factor 7 (IRF7) and a significant lower abundance of IRF7 was found in patients with COVID-19 compared to influenza (Supplementary figure 5). The second differentially activated TopMD pathway had peak gene NDUFAB1; named by GO analysis as mitochondrial complex I assembly model OXPHOS system WP4324 (adjusted p-value 2.81×10⁻⁶⁶). The third differentially activated TopMD pathway was named by GO analysis as proteasome degradation WP183 (adjusted p-value, 1.46×10⁻⁶⁴), with PSMD14 as the peak gene (Figure 5 with full detail in Additional file 3 and the global map of differentially activated pathways available online (65)).

• Download figure
• Open in new tab

Figure 5: Differentially activated pathways between hospitalised patients with COVID-19 or influenza identified with topological analysis.

The difference (Log2 fold change) in patients with COVID-19 compared to patients with influenza is plotted for the top 20 genes of the 1^st, 2^nd and 3^rd TopMD pathways.

Cell subsets supporting innate and adaptive immune response differences

Analysis of the blood transcriptome can be used to predict the immune cells present (64). Levels of different predicted cell types were assessed to determine whether there were differences in immune system associated cells between patients with COVID-19 or influenza (Figure 6). Statistical testing was done on cell type levels identified with CIBERSORTx. M0 macrophages (p-value 3.63×10⁻⁰⁶), plasma cells (p-value 5.05×10⁻⁰⁴), cytotoxic CD8+ T cells (p-value 4.58×10⁻⁰³), regulatory T cells (p-value 7.30×10⁻⁰³) and resting natural killer cell (p-value 8.90×10⁻⁰³) were found to be significantly higher in COVID-19 patients, while in influenza patients activated dendritic cells (p-value 2.23×10⁻⁰²) were significantly higher.

• Download figure
• Open in new tab

Figure 6: Increase of predicted cells associated with an innate and adaptive immune response in COVID-19 patients.

M0 macrophages, plasma cells, cytotoxic CD8+ T cells, regulatory T cells and resting natural killer (NK) cells were found to be significantly higher in COVID-19 patients. In influenza patients a significantly higher proportion of activated dendritic cells was detected.

Predicted cell type levels between COVID-19 survivors and non-survivors indicated an increase of neutrophils (p-value 2.84×10⁻⁰⁴) in patients who died of COVID-19. In contrast, an increase of naïve CD4+ T cells (p-value 1.92×10⁻⁰³), M0 macrophages (p-value 1.20×10⁻⁰²), M2 macrophages (p-value 1.48×10⁻⁰²), naïve B cells (p-value 1.57×10⁻⁰²) and naïve cytotoxic CD8+ T cells (p-value 2.31×10⁻⁰²), were identified in patients who went on to survive COVID-19 (Figure 7).

• Download figure
• Open in new tab

Figure 7: Differences in immune response indicated by predicted cell types in COVID-19 survivors and non-survivors.

A) A statistically significant higher count of neutrophils in COVID-19 patients who died after 30 days indicating the presence of an elevated innate immune response. B) Adaptive immune response in COVID-19 survivors as can be seen by the statistically significant higher count of naïve B cells, and CD4+ and CD8+ T cells.

Efficient adaptive immune response associates with COVID-19 survival

As already noted, a high abundance of the GO biological process ‘adaptive immune response’ related transcripts, mostly immunoglobulin genes, was associated with COVID-19 survival (Figure 4). Here a direct assessment was done of the blood transcriptome differences between patients who, at 30 days after hospital admission, survived or who died of COVID-19. A total of 23,850 abundance measures of gene transcripts were obtained after filtering out transcripts with low counts, of which 6.645 transcripts were found to be significant (FDR p-value < 0.05) of which, 537 transcripts exceeded a log2 fold change of < −1 or > 1, with 265 transcripts showing higher abundance in patients who survived COVID-19 and 272 transcripts showing a higher abundance in patients who died of COVID-19 (Figure 8A and Additional file 4).

• Download figure
• Open in new tab

Figure 8: Increased innate immune response in COVID-19 non-survivors and increased adaptive immune response in COVID-19 survivors.

A) Volcano plot of transcripts falling within the threshold values (FDR < 0.05 and log fold change < −1 or >1) which were used for enrichment analysis with ToppGene. B) Enrichment analysis identified an increased innate immune response in patients who died of COVID-19 after 30 days of hospital admission. Percentage in annotation is the ratio of the input query genes overlapping with the genes in the pathway database. C) An enrichment of adaptive immune response related pathways was detected in patients with COVID-19 who were still alive 30 days after hospital admission. Percentage in annotation is the ratio of the input query genes overlapping with the genes in the pathway database. D) Increase of T cell and B cell proliferation in COVID-19 survivors (paired non-parametric T-test).

In patients who died of COVID-19 an enrichment for biological processes involved with an inflammatory response including interleukin signalling and neutrophil activation and degranulation was detected (Figure 8B). While in COVID-19 survivors biological processes involved with the adaptive immune system including complement activation, B cell mediated immunity and a humoral immune response mediated by circulating immunoglobulins was found to be enriched (Figure 8C). Additionally, transcript abundances associated with T cell and B cell proliferation were significantly higher in COVID-19 survivors (p-value < 1.0×10⁻⁰⁴, paired non-parametric T-test) (Figure 8D).

Immune signatures as predictors of COVID-19 outcome

Distinct immune signature genes were selected and assessed for their prediction accuracy in stratifying patients with COVID-19 for disease outcome, fatality or survival. A signature of 47 genes was identified (Figure 9A), representative of the four biggest clusters of genes associated with either patients with COVID-19 who survived or died. The associated GO biological process terms were ‘humoral immune response mediated by circulating immunoglobulin’ (FDR p-value 2.23×10⁻⁴⁶), ‘nucleosome assembly’ (FDR p-value 5.46×10⁻¹⁹), ‘regulation of T-helper 1 cell cytokine production’ (FDR p-value 4.24×10⁻⁰³) and ‘regulation of T cell activation’ (FDR p-value 4.51×10⁻⁰⁴) (Supplementary figure 6). This was highly predictive for outcome, with a maximum specificity of 75% and sensitivity of 93% (Figure 9B and Table 4).

• Download figure
• Open in new tab

Figure 9: Receiver Operating Characteristic (ROC) curves showing prediction accuracy COVID-19 survivors and non-survivors.

A) Genes identified with EdgeR and gene co-expression analysis and used for subsequent modelling. B) ROC curves according to the three models used (Boosted Logistic Regression (LogitBoost), Bayesian Generalised Linear (Bayesglm) and RandomForest (rf)).