In vivo antiviral host transcriptional response to SARS-CoV-2 by viral load, sex, and age

Abstract

Despite limited genomic diversity, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has shown a wide range of clinical manifestations in different patient populations. The mechanisms behind these host differences are still unclear. Here, we examined host response gene expression across infection status, viral load, age, and sex among shotgun RNA sequencing profiles of nasopharyngeal (NP) swabs from 430 individuals with PCR-confirmed SARS-CoV-2 and 54 negative controls. SARS-CoV-2 induced a strong antiviral response with up-regulation of antiviral factors such as OAS1-3 and IFIT1-3 and T helper type 1 (Th1) chemokines CXCL9/10/11, as well as a reduction in transcription of ribosomal proteins. SARS-CoV-2 culture in human airway epithelial (HAE) cultures replicated the in vivo antiviral host response 7 days post infection, with no induction of interferon-stimulated genes after 3 days. Patient-matched longitudinal specimens (mean elapsed time = 6.3 days) demonstrated reduction in interferon-induced transcription, recovery of transcription of ribosomal proteins, and initiation of wound healing and humoral immune responses. Expression of interferon-responsive genes, including ACE2, increased as a function of viral load, while transcripts for B cell-specific proteins and neutrophil chemokines were elevated in patients with lower viral load. Older individuals had reduced expression of the Th1 chemokines CXCL9/10/11 and their cognate receptor CXCR3, as well as CD8A and granzyme B, suggesting deficiencies in trafficking and/or function of cytotoxic T cells and natural killer (NK) cells. Relative to females, males had reduced B cell-specific and NK cell-specific transcripts and an increase in inhibitors of nuclear factor kappa-B (NF-κB) signaling, possibly inappropriately throttling antiviral responses. Collectively, our data demonstrate that host responses to SARS-CoV-2 are dependent on viral load and infection time course, with observed differences due to age and sex that may contribute to disease severity.

Fig 1. DE genes in SARS-CoV-2 NP swabs.

(A) Clustering of samples based on 50 genes with the lowest adjusted p-value. Log2 fold changes relative to gene mean are displayed by color. (B) Volcano plot of 15 most up-regulated and 15 most down-regulated genes in SARS-CoV-2 positive samples relative to negative by log2 fold change. Red color indicates genes with log2 fold change > |1.5| and adjusted p < 0.05. (C). Significant (FDR < 0.05) pathways affected by SARS-CoV-2 infection identified by GSEA. Raw data available in the GEO Repository, accession GSE152075. DE, differentially expressed; FDR, false discovery rate; GEO, Gene Expression Omnibus; GSEA, Gene Set Enrichment Analysis; NP, nasopharyngeal; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

Fig 2. Differences in gene expression by SARS-CoV-2 viral load.

(A) Violin plots of select genes by viral load. Statistical significance between low and high viral load calculated by Mann Whitney U test, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. (B) Volcano plot of 15 most up-regulated and 15 most down-regulated genes in SARS-CoV-2 high viral load samples relative to low viral load by log2 fold change. Red color indicates genes with log2 fold change > |1.5| and adjusted p < 0.05. (C) Proportion of cell types as a total of all immune cells, by CIBERSORTx. Significant differences in proportion of each cell type determined by t test, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. (D) Violin plots of B cell transcripts and neutrophil chemokine transcripts by viral load. Statistical significance between low and high viral load calculated by Mann Whitney U test, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. Raw data available in the GEO Repository, accession GSE152075. Ct, cycle threshold; GEO, Gene Expression Omnibus; N1, SARS-CoV-2 nucleocapsid gene region 1; NK, natural killer; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2; Treg, regulatory T cell.

Fig 3. Consensus genes induced upon SARS-CoV-2 expression.

(A) Venn diagram of DE genes in SARS-CoV-2-positive versus -negative, high versus low viral load, and top 100 genes with the highest absolute log2 fold change in infected versus uninfected HAE. Consensus set of 19 genes (SuperExact Test, p = 4.81 × 10⁻⁶⁹) DE in all 3 analyses are shown, with log2 fold change values relative to uninfected HAE (for day 3 and day 7 post infection), SARS-CoV-2-negative NP swabs (for SARS-CoV-2-positive NP swabs), or low viral load (for high SARS-CoV-2 viral load samples). SARS-CoV-2 reads at day 3 and 7 post infection were 0.3% and 5.3%, respectively. (B) Top 20 DisGeNET terms for which SARS-CoV-2 cell-intrinsic antiviral response consensus genes are overrepresented. “Number Enriched” is the number of SARS-CoV-2 consensus genes that belong to each disease term. (C) Interaction network of SARS-CoV-2 consensus genes for top 5 most similar diseases identified in panel B. Size of disease node represents the number of genes enriched, and fold change is the log2 fold change seen in SARS-CoV-2-positive versus -negative NP swabs. Raw data available in the GEO Repository, accession GSE154768 (HAE samples) or GSE152075 (NP swabs). DE, differentially expressed; GEO, Gene Expression Omnibus; HAE, human airway epithelial cells; NP, nasopharyngeal; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

Fig 4. DE genes in patient-matched longitudinal samples.

(A) Patient demographics information for longitudinal samples. (B) Top 20 Biological Process GO terms for which longitudinal DE genes are overrepresented. “Number Enriched” is the number of DE genes that belong to each GO term. (C) Log2 fold changes for DE genes in Humoral Immune Response and Wound Healing GO terms, consensus antiviral SARS-CoV-2 genes, and ribosomal proteins. Grey bars: p_adj < 0.1, white bars: p_adj > 0.1. Raw data available in the GEO Repository, accession GSE154769. Ct, cycle threshold; DE, differentially expressed; GEO, Gene Expression Omnibus; GO, Gene Ontology; N1, SARS-CoV-2 nucleocapsid gene region 1; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.

Fig 5. Age and sex cause differences in gene expression upon SARS-CoV-2 infection.

(A) N1 Ct values by age group. No significant differences between were observed by Kruskal-Wallis ANOVA. (B) N1 Ct values by sex. No significant difference between groups was observed by t test. (C) Gene expression differences by age and viral load. Significance by Mann Whitney U test between SARS-CoV-2-positive samples aged >60 and <60 is shown, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. (D) Sex-modulated DE genes (p_adj < 0.1) upon SARS-CoV-2 infection. Genes elevated in females are shown as negative log2 fold changes, and those elevated in males as positive log2 fold changes. Raw data available in the GEO Repository, accession GSE152075. Ct, cycle threshold; DE, differentially expressed; GEO, Gene Expression Omnibus; N1, SARS-CoV-2 nucleocapsid gene region 1; NK, natural killer; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2.