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

  1. Matias A Medina
  2. Francisco Fuentes-Villalobos
  3. Claudio Quevedo
  4. Felipe Aguilera
  5. Raul Riquelme
  6. Maria Luisa Rioseco
  7. Sebastian Barria
  8. Yazmin Pinos
  9. Mario Calvo
  10. Ian Burbulis
  11. Camila Kossack
  12. Raymond A Alvarez
  13. Jose Luis Garrido  Is a corresponding author
  14. Maria Ines Barria  Is a corresponding author
  1. Facultad de Medicina y Ciencia, Universidad San Sebastián, Chile
  2. Departamento de Microbiología, Facultad de Ciencias Biológicas, Universidad de Concepción, Chile
  3. Departamento de Bioquímica y Biología Molecular, Facultad de Ciencias Biológicas, Universidad de Concepción, Chile
  4. Hospital Dr. Eduardo Schütz Schroeder, Chile
  5. Hospital Base San José, Chile
  6. Instituto de Medicina, Facultad de Medicina, Universidad Austral, Chile
  7. Division of Infectious Diseases, Department of Medicine, Immunology Institute, Icahn School of Medicine at Mount Sinai, United States
5 figures, 1 table and 1 additional file

Figures

Figure 1 with 2 supplements
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.

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.

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 3—source data 1

Enriched pathways from genes upregulated in severe COVID-19 patients by longitudinal analysis for Figure 3A, C, E.

https://cdn.elifesciences.org/articles/94242/elife-94242-fig3-data1-v1.docx
Figure 3—source data 2

Enriched pathways from genes upregulated in severe COVID-19 patients by longitudinal analysis for Figure 3B, D, F.

https://cdn.elifesciences.org/articles/94242/elife-94242-fig3-data2-v1.docx
Figure 4 with 4 supplements
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—source data 1

Kyoto encyclopedia of gene and genomes (KEGG) graph shows genes with differential expression found in the pairwise comparison (D0 vs D7) from natural killer cell-mediated cytotoxicity pathway for upregulated genes in mild COVID-19 patients at D0 (A) and D7 (B), from Th1 and Th2 cell differentiation pathway at D0 (C) and D7 (D), and from cytokine–cytokine receptor interaction pathway at D0 (E).

Red boxes depict upregulated genes, whereas green boxes depict genes without significant differential gene expression within each KEGG pathway.

https://cdn.elifesciences.org/articles/94242/elife-94242-fig4-data1-v1.pdf
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 with 2 supplements
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.

Tables

Table 1
Clinical characteristics of the cohort.
Clinical characteristicsMildSevereAll
Sex, n (female/male)(2/2)(2/2)(4/4)
Total, n448
Median age, years ± SD39.0 ± 3.946.7 ± 842.9 ± 7
Days from onset of symptoms to recruitment, median ± SD1.2 ± 1.310.0 ± 1.85 ± 4.7
Days from COVID-19 diagnosis to recruitment, median ± SD3.0 ± 0.78.5 ± 2.16 ± 3.4
Symptoms, n (%)
Fever1 (25%)3 (75%)4 (50%)
Chills1 (25%)2 (50%)3 (37.5%)
Fever feeling1 (25%)2 (50%)3 (37.5%)
Odynophagia2 (50%)1 (25%)3 (37.5%)
Cough1 (25%)4 (100%)5 (62.5%)
Expectoration---
Dyspnea-4 (100%)4 (50%)
Thoraxic pain1 (25%)1 (25%)2 (25%)
Diarrhea1 (25%)2 (50%)3 (37.5%)
Anosmia1 (25%)-1 (12.5%)
Ageusia1 (25%)-1 (12.5%)
Myalgia1 (25%)2 (50%)3 (37.5%)
Headache2 (50%)1 (25%)3 (37.5%)
Treatment
Hospitalization, n (%)-4 (100%)4 (50%)
Intubation, n (%)-4 (100%)4 (50%)
Mechanical ventilation, n (%)-4 (100%)4 (50%)
Samples, n121022

Additional files

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Matias A Medina
  2. Francisco Fuentes-Villalobos
  3. Claudio Quevedo
  4. Felipe Aguilera
  5. Raul Riquelme
  6. Maria Luisa Rioseco
  7. Sebastian Barria
  8. Yazmin Pinos
  9. Mario Calvo
  10. Ian Burbulis
  11. Camila Kossack
  12. Raymond A Alvarez
  13. Jose Luis Garrido
  14. Maria Ines Barria
(2024)
Longitudinal transcriptional changes reveal genes from the natural killer cell-mediated cytotoxicity pathway as critical players underlying COVID-19 progression
eLife 13:RP94242.
https://doi.org/10.7554/eLife.94242.3