Longitudinal host transcriptional responses to SARS-CoV-2 infection in adults with extremely high viral load

Abstract

Current understanding of viral dynamics of SARS-CoV-2 and host responses driving the pathogenic mechanisms in COVID-19 is rapidly evolving. Here, we conducted a longitudinal study to investigate gene expression patterns during acute SARS-CoV-2 illness. Cases included SARS-CoV-2 infected individuals with extremely high viral loads early in their illness, individuals having low SARS-CoV-2 viral loads early in their infection, and individuals testing negative for SARS-CoV-2. We could identify widespread transcriptional host responses to SARS-CoV-2 infection that were initially most strongly manifested in patients with extremely high initial viral loads, then attenuating within the patient over time as viral loads decreased. Genes correlated with SARS-CoV-2 viral load over time were similarly differentially expressed across independent datasets of SARS-CoV-2 infected lung and upper airway cells, from both in vitro systems and patient samples. We also generated expression data on the human nose organoid model during SARS-CoV-2 infection. The human nose organoid-generated host transcriptional response captured many aspects of responses observed in the above patient samples, while suggesting the existence of distinct host responses to SARS-CoV-2 depending on the cellular context, involving both epithelial and cellular immune responses. Our findings provide a catalog of SARS-CoV-2 host response genes changing over time and magnitude of these host responses were significantly correlated to viral load.

Fig 1. Differential gene sets associated with the transcriptional host response to SARS-CoV-2 infection across serially collected samples.

(a) For each of 20000 genes, expression (log2 FPKM) was correlated with viral load (inverse correlation with Ct value) across 44 samples from 20 subjects. Numbers of statistically significant genes (by Pearson’s) at both p<0.01 and p<0.001 significance levels are represented, as compared to the chance expected by multiple testing [25]. (b) Numbers of differential genes (p<0.01, t-test) when comparing: 1) Visit 1 samples from the extremely high viral load group (n = 8 samples from eight subjects) with the samples in the negative group (n = 4); samples at the latest time points for each of the subjects from the extremely high viral load group (n = 8 samples) with the samples in the negative group; samples from the low viral load group (n = 8 samples from eight subjects, using earliest time point) with the samples in the negative group. Chance expected genes at p<0.01 due to multiple testing would be on the order of 200 [25]. (c) Heat map comparing differential patterns across the three comparisons from part b, for the 1357 genes significant (p<0.01) for any comparison. Row along the bottom indicated which genes were positively correlated with viral load (p<0.01) across all 44 samples (from part a), and which genes have Gene Ontology (GO) annotation [26] ‘response to virus’. Visualization using heat maps was performed using JavaTreeview (version 1.1.6r4) [27].

Fig 2. UpSet plot summarising key differentially expressed gene trends and Gene ontology analysis of differentially expressed genes present only in extremely high viral load group.

(a) The plot depicts common and unique genes shared between extremely high viral load group and low viral load group as compared to SARS-CoV-2 negative group groups. Set size indicates number of differentially expressed genes in each comparison. Intersection size is the number of statistically significant (FDR < 0.05) differentially expressed genes in designated sets or groups. Connected circles at the bottom of the plot indicate an intersection of differentially expressed genes between groups. (b)The top10 enriched GO terms for biological processes altered in extremely high viral load group in shown. Significantly enriched GO terms with a minimum three enriched genes were ranked by significance. The x-axis denoting the negative log fold change of significance. Dotted red line depicts the significance threshold of FDR <0.05.

Fig 3. Differential expression patterns and functional gene groups associated with SARS-CoV-2 viral load across serially collected samples.

(a) Across 44 MT swab samples representing 20 subjects, differential gene expression patterns for the set of 112 genes significantly correlated with SARS-CoV-2 viral load (i.e., inversely correlated with Ct value) at p<0.001 (Pearson’s) are represented. Heat map contrast (bright yellow/blue) is 3-fold change from the average of the samples from the low viral load group. Genes listed off to the right have GO annotation ‘response to virus’. Extremely high viral load, Ct<20. Visualization using heat maps was performed using JavaTreeview (version 1.1.6r4) [27] (b) Selected significantly enriched GO terms [26] within the genes over-expressed with SARS-CoV-2 viral load (p<0.01, Pearson’s). For each GO term, enrichment p-values and numbers of genes in the SARS-CoV-2-associated gene set are indicated. Enrichment p-values by one-sided Fisher’s exact test.

Fig 4. Genes correlated with SARS-CoV-2 viral load over time are similarly expressed in independent datasets of SARS-CoV-2 infected lung and upper airway cells.

(a) Differential expression patterns for the 112 genes correlated with SARS-CoV-2 viral load across our serial sampling cohort (p<0.001, from Fig 3A) were examined in two independent RNA-seq datasets of SARS-CoV-2 infection: one of lung cancer cell lines (A549 and Calu-3) infected with SARS-CoV-2 at multiplicity-of-Infection (MOI) of 2 for 24 hours (in biological triplicate)[28], and one of nasopharyngeal/oropharyngeal samples in 238 patients with COVID-19, other viral, or non-viral acute respiratory illnesses [29]. Gene order is the same across all datasets. Heat map contrast (bright yellow/blue) is 3-fold change from the corresponding comparison group (serial sampling dataset, average of the samples from the low viral load group; lung cancer cell line dataset, average of corresponding mock control group; Mick et al. dataset, average of “no virus” samples). Visualization using heat maps was performed using JavaTreeview (version 1.1.6r4) [27] (b) Venn diagram representing the gene set overlaps among the genes increased with SARS-CoV-2 infection in each of the three RNA-seq datasets from part a (with Calu-3 lung cancer cell line being considered here over A549). A p-value cutoff of p<0.01 was used to define top genes for each dataset (serial MT swab and Mick et al. nasopharyngeal/oropharyngeal datasets, Pearson’s correlation with viral load; Calu-3 dataset, t-test). Gene set enrichment p-values by one-sided Fisher’s exact test. Genes overlapping between all three datasets are listed. (c) Similar to part b, but for genes decreased with SARS-CoV-2 infection.

Fig 5. Differential expression patterns and functional gene groups associated with SARS-CoV-2 infection of nose organoids.

(a) HNO204 human nose organoids were infected with SARS-CoV-2 at an MOI of 0.01, and samples at 6hrs, 72hrs, and 6 days post infection were profiled for gene expression. Differential expression patterns for the top 867 genes over-expressed in HNO204 with SARS-CoV-2 infection (p<0.05, t-test) are represented here. Next to the HNO204 dataset are the corresponding patterns for independent RNA-seq datasets of SARS-CoV-2 infection: lung cancer cell lines (A549 and Calu-3)²⁵, our serially collected MT swab samples from patients, and nasopharyngeal/oropharyngeal samples from Mick et al [29]. Gene order is the same across all datasets. Heat map contrast (bright yellow/blue) is 3-fold change from the corresponding comparison group. Regarding the Mick et al. dataset, the ordering of the samples is the same as that for Fig 4A, and the SARS-CoV-2 viral load plot represents log2 RPM, with values ranging from 0 to 19.6. Visualization using heat maps was performed using JavaTreeview (version 1.1.6r4) [27] (b) Venn diagram representing the gene set overlaps among the genes increased with SARS-CoV-2 infection in each of the following RNA-seq datasets: HNO204, serial MT swab, and Calu-3 lung cancer cell line. Gene set enrichment p-values by one-sided Fisher’s exact test. Genes overlapping between HNO204 and serial MT swab datasets are listed. (c) Similar to part b, but for genes decreased with SARS-CoV-2 infection. (d) Selected significantly enriched GO terms ⁴¹ within the genes over-expressed with SARS-CoV-2 infection in HNO204 (p<0.05, t-test). For each GO term, enrichment p-values and numbers of genes in the SARS-CoV-2-associated gene set are indicated. Enrichment p-values by one-sided Fisher’s exact test.