Longitudinal transcriptional changes reveal genes from the natural killer cell-mediated cytotoxicity pathway as critical players underlying COVID-19 progression

Abstract

Patients present a wide range of clinical severities in response severe acute respiratory syndrome coronavirus 2 infection, but the underlying molecular and cellular reasons why clinical outcomes vary so greatly within the population remains unknown. Here, we report that negative clinical outcomes in severely ill patients were associated with divergent RNA transcriptome profiles in peripheral immune cells compared with mild cases during the first weeks after disease onset. Protein-protein interaction analysis indicated that early-responding cytotoxic natural killer cells were associated with an effective clearance of the virus and a less severe outcome. This innate immune response was associated with the activation of select cytokine-cytokine receptor pathways and robust Th1/Th2 cell differentiation profiles. In contrast, severely ill patients exhibited a dysregulation between innate and adaptive responses affiliated with divergent Th1/Th2 profiles and negative outcomes. This knowledge forms the basis of clinical triage that may be used to preemptively detect high-risk patients before life-threatening outcomes ensue.

Keywords: COVID-19; NK; PBMC RNA-seq; human; infectious disease; microbiology.

Figure 1.. Clinical profile and gene expression patterns for mild and severe COVID-19 patients during 28 days.

(A) Longitudinal sampling schedule for severe and mild COVID-19 donors of peripheral blood (22 samples in total from 8 donors). Three sampling times (black dots) are displayed with respect to the recruitment day (D0). In addition, panel A shows the diagnosis day with positive PCR (green star), symptoms onset day (orange star), and hospitalization day (red cross). (B) The COVID-19 progression according to the WHO ordinal scale describes the temporal disease severity for each donor. Severe patients (S1–S4) are displayed in lines with colors scaling from green to red, while mild patients (M1–M4) are shown in blue line colors. The X-axis represents the days relative to the recruitment day. The vertical dot lines indicate the three sampling times used in this study (days 0, 7, and 28). (C) Principal component analysis plot based on gene expression profiles for mild (blue) and severe (red) COVID-19 patients and grouped by their sampling time (0, 7, and 28 days after recruitment). In addition, peripheral blood samples from two healthy donors are shown in green. PC = principal component. (D) Heatmap of the 100 most significant differentially expressed genes related to the COVID-19 disease progression. At the bottom, each column corresponds to the sampling points (0, 7, and 28 days since recruitment) of mild and severe patients. Genes are displayed as horizontal rows and are clustered by the similarity of expression profiles, represented by the dendrogram to the left of the heatmap. Red indicates higher expression, while blue means lower expression represented by the z-score of normalized read counts. Some COVID-19 severity-associated genes previously reported are indicated to the right of the heatmap in red.

Figure 1—figure supplement 1.. Heatmap of temporally and differentially expressed genes over the course of COVID-19 progression.

At the top, each column corresponds to the sampling points (0, 7, and 28 days since recruitment) of mild and severe patients. Genes are displayed as horizontal rows and are clustered by the similarity of expression profiles, represented by the dendrogram to the left of the heatmap. To the right of the heatmap, yellowish color indicates higher expression, while bluish color means lower expression represented by the z-score of normalized read counts.

Figure 1—figure supplement 2.. Volcano plot depicting pairwise gene expression comparisons for detected differentially expressed genes (DEGs) between mild and severe COVID-19 patients at D0, D7, and D28 (A, B, and C, respectively).

Red or blue indicates genes that were significantly over- or sub-expressed at a particular sampling time, based on filtering by FDR ≤0.05 and the absolute value of −log10 (FC) (≥2.0). The remaining genes that do not show differential expression are indicated in gray. FDR = false discovery rate; FC = fold change.

Figure 2.. Gene set enrichment analysis (GSEA) shows biological processes and enriched pathways associated to innate and adaptive immune response after severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection.

(A, B) Bubble plots show the biological processes (BP) of enriched genes for differentially expressed genes (DEGs) between mild and severe for D0 and D7 after recruitment, respectively. (C, D) Bubble plots show the Kyoto encyclopedia of gene and genomes (KEGG) enriched pathways of enriched genes for DEGs between mild and severe for D0 and D7, respectively. The Count (black circles) represents the number of genes included on each set. Generatio is the ratio between number of genes found in the set and total genes of set. The scale-color bar indicates the p-adjustment of each BP or KEGG pathway.

Figure 3.. Biological processes and Kyoto encyclopedia of gene and genomes (KEGG) pathway for genes with differential expression levels over time in mild versus severe patients.

Bar plots show gene ontology analysis for biological processes in light-blue/dark-blue and KEGG pathway in light-green/dark-green, using upregulated genes (A, C, E) and downregulated genes (B, D, F). The size of each bar is according to its −log10 p-value, and name pathways with asterisk indicates a q-value ≤0.05.

Figure 4.. Gene function network for natural killer (NK) cell hub-genes with differential expression levels between mild and severe patients during acute phase of COVID-19.

(A) Kyoto encyclopedia of gene and genomes (KEGG) pathway of NK cell-mediated cytotoxicity represents the set NK cell hub-genes upregulated (red-boxes) in mild versus severe patients. The green boxes correspond to genes without differential expression. (B) Heatmap shows the differential expression levels of NK cell hub-genes over time (D0, D7, and D28 after recruitment) separated by mild and severe groups. The expression levels are represented by the z-score of normalized counts. Dendrogram shows the hierarchical clustering of genes. (C) Protein–protein interaction (PPI) network for upregulated genes during the acute phase in mild patients. The network corresponds to the principal clusters with more interaction between proteins and highlights the three most represented pathways: Cytokine–cytokine receptor interaction (blue); NK cell-mediated cytotoxicity (red); and Th1 and Th2 cell differentiation (yellow). (D) Time course expression levels for the main protein nodes identified in PPI network during the acute phase of COVID-19. The trajectories of these genes are graphed as days after recruitment (D0, D7, and D28) for mild (red triangle) and severe (green circle) groups and their enrichment is represented by the z-score of normalized counts.

Figure 4—figure supplement 1.. Kyoto encyclopedia of gene and genomes (KEGG) graphs show genes differentially expressed over time of the Th1 and Th2 cell differentiation pathway (A) and the cytokine–cytokine receptor interaction pathway (B).

Red boxes depict upregulated genes in mild COVID-19 patients during the acute phase of the disease, whereas green boxes depict genes without significant differential gene expression within the KEGG pathways.

Figure 4—figure supplement 2.. Gene ontology (GO) and Kyoto encyclopedia of gene and genomes (KEGG) and REACTOME pathway analyses of differentially expressed genes (DEGs) found in the pairwise comparison between D0 and D7 of COVID-19 infection.

Bar graphs showing the enrichment of GO biological processes (in blue color), GO molecular functions (in purple color), KEGG pathways (in dark green color), and Reactome pathways (in light green color) between mild and severe COVID-19 patients. Enrichment results are sorted by −log10(p-value) (higher on top) with a cut-off for DEGs ≥4 fold-change considering the upregulated genes at D0 (A), upregulated genes at D7 (B), downregulated genes at D0 (C), and downregulated genes at D7 (D).

Figure 4—figure supplement 3.. Protein–protein interaction (PPI) network graphs show the upregulated genes found in the pairwise comparison (D0 and D7) in mild versus severe COVID-19 patients.

(A) Topological representation of the PPI network of upregulated genes in mild patients on D0. (B) Topological representation of the PPI network of upregulated genes in mild patients on D7. Some nodes are color-coded to highlight proteins involved in the following pathways: Cytokine–cytokine receptor interaction (blue), natural killer (NK) cell-mediated cytotoxicity (red); and Th1 and Th2 cell differentiation (yellow).

Figure 4—figure supplement 4.. Gene ontology (GO) and networks of enriched pathway analyses of differentially expressed genes (DEGs) found in the pairwise comparison between D7 of mild patients and D0 of severe patients with COVID-19.

Bar graphs show the enrichment of GO terms for biological processes while the networks reveal Kyoto encyclopedia of gene and genomes (KEGG) enriched pathways for upregulated genes (A, B) and downregulated genes (C, D) in mild patients, respectively. The size of each bar is according to its −log10 p-value, and name GO terms with asterisk indicates a q-value ≤0.05.

Figure 5.. Gene co-expression network analysis among the longitudinal transcriptomic profiling.

(A) Gene hierarchical clustering dendrogram of 10 detected modules based on Topological Overlap Matrices (TOM) measure. The branches and color bands represent the assigned module. (B) Module–trait relationships (MTRs) between detected modules and disease severity of COVID-19. MTRs are obtained by calculating Pearson correlation between the traits and the module eigengenes. The red and blue colors indicate strong positive or negative correlations, respectively. Rows represent module eigengene (ME) and columns indicate the disease severity of COVID-19, where correlations from each patient are shown from left to right according to their sampling time (D0, D7, and D28). (C) Gene ontology (GO) enrichment analysis of genes in the blue, brown, and turquoise modules. The color in the bar graphs refers to the module eigengene (ME). Enrichment results are sorted by −log10(p-value) (higher on top) of each biological process GO term. (D–F) Network visualization of enriched pathways based on co-expression genes from blue, brown, and turquoise modules, respectively.

Figure 5—figure supplement 1.. Gene ontology (GO) analysis of the seven smallest modules of co-expression.

The genes from the modules red, black, magenta, green, pink, yellow, and gray were used to study GO biological process terms. GO terms with statistically significant q-values are marked with an asterisk.

Figure 5—figure supplement 2.. Network for genes of enriched pathways from the seven smallest modules of co-expression.

Genes from modules red (A), black (B), pink (C), yellow (D), green (E), magenta (F), and gray (G) were used for analysis and the networks show the most enriched pathways.