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.

About the Authors:

Nicole A. P. Lieberman

Roles Conceptualization, Formal analysis, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0001-9334-8150

Vikas Peddu

Roles Data curation, Formal analysis

Affiliation: Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-0366-3956

Hong Xie

Roles Methodology

Affiliation: Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America

Lasata Shrestha

Roles Methodology

Affiliation: Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-9760-7348

Meei-Li Huang

Roles Methodology, Supervision

Affiliation: Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America

Megan C. Mears

Roles Methodology

Affiliations Galveston National Laboratory, University of Texas Medical Branch, Galveston, Texas, United States of America, Department of Experimental Pathology, University of Texas Medical Branch, Galveston, Texas, United States of America

ORCID logo http://orcid.org/0000-0003-1293-7696

Maria N. Cajimat

Roles Methodology

Affiliations Galveston National Laboratory, University of Texas Medical Branch, Galveston, Texas, United States of America, Department of Experimental Pathology, University of Texas Medical Branch, Galveston, Texas, United States of America

Dennis A. Bente

Roles Supervision

Affiliations Galveston National Laboratory, University of Texas Medical Branch, Galveston, Texas, United States of America, Department of Microbiology and Immunology, University of Texas Medical Branch, Galveston, Texas, United States of America

ORCID logo http://orcid.org/0000-0002-5150-7368

Pei-Yong Shi

Roles Resources

Affiliations Galveston National Laboratory, University of Texas Medical Branch, Galveston, Texas, United States of America, Department of Biochemistry and Molecular Biology, University of Texas Medical Branch, Galveston, Texas, United States of America

Francesca Bovier

Roles Data curation, Methodology

Affiliation: Center for Host–Pathogen Interaction, Columbia University Medical Center, New York, New York, United States of America

Pavitra Roychoudhury

Roles Data curation

Affiliations Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America, Vaccine and Infectious Disease Division, Fred Hutchinson Cancer Research Center, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-4567-8232

Keith R. Jerome

Roles Funding acquisition, Writing – review & editing

Affiliations Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America, Vaccine and Infectious Disease Division, Fred Hutchinson Cancer Research Center, Seattle, Washington, United States of America

Anne Moscona

Roles Supervision

Affiliations Center for Host–Pathogen Interaction, Columbia University Medical Center, New York, New York, United States of America, Department of Pediatrics, Columbia University Medical Center, New York, New York, United States of America, Department of Microbiology & Immunology, Columbia University Medical Center, New York, New York, United States of America, Department of Physiology & Cellular Biophysics, Columbia University Medical Center, New York, New York, United States of America

Matteo Porotto

Roles Funding acquisition, Supervision

Affiliations Center for Host–Pathogen Interaction, Columbia University Medical Center, New York, New York, United States of America, Department of Pediatrics, Columbia University Medical Center, New York, New York, United States of America, Department of Experimental Medicine, University of Campania “Luigi Vanvitelli,” Caserta, Italy

Alexander L. Greninger

Roles Conceptualization, Funding acquisition, Supervision, Writing – review & editing

* E-mail: [email protected]

Affiliations Department of Laboratory Medicine and Pathology, University of Washington School of Medicine, Seattle, Washington, United States of America, Vaccine and Infectious Disease Division, Fred Hutchinson Cancer Research Center, Seattle, Washington, United States of America

ORCID logo http://orcid.org/0000-0002-7443-0527

Introduction

The novel coronavirus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) that emerged in late 2019 from Wuhan, China, has rapidly spread throughout the world, causing more than 6 million cases and 400,000 deaths globally as of June 2020. Coronavirus disease 2019 (COVID-19) morbidity and mortality has been overwhelmingly concentrated in elderly individuals and those with preexisting comorbidities [1]. In older individuals, immunosenescence and dysregulated antiviral responses due to viral chronic low-grade age-related inflammation may play an important role [2], as has been proposed for influenza [3]. Males are known to be generally more susceptible to infectious disease than females [4], and severe acute respiratory syndrome coronavirus (SARS-CoV)-infected male mice had increased infiltration of inflammatory macrophages into their lungs, leading to a deleterious inflammatory response [5]. Accordingly, systemic inflammatory markers such as neutrophil-to-lymphocyte ratio and C-reactive protein were elevated in men who died of SARS-CoV-2 [6]. However, the mechanisms behind increased mortality among older adults and males with COVID-19 remain speculative.

Entry of SARS-CoV-2 into host cells depends on binding to the receptor angiotensin-converting enzyme 2 (ACE2) [7], expressed at a high level in the nasal epithelium [8], then further induced upon exposure to interferon [9], suggesting a mechanism by which SARS-CoV-2 exploits host antiviral responses. SARS-CoV antagonizes initial viral detection and interferon responses by an as-yet unknown mechanism [10,11]. SARS-CoV-2 may employ similar mechanisms, as infections of bronchial epithelial cells at low multiplicity of infection (MOI) do not result in extensive transcription of interferon-stimulated genes (ISGs) at 24 hours post infection [12]. An important consequence of these observations is that SARS-CoV-2 viral load and transmissibility peaks at the time of symptom onset [13,14]. The temporal relationship between viral load and host gene expression has not been fully explored.

In the United States, diagnostic testing is generally performed on nasopharyngeal (NP) swabs, from which SARS-CoV-2 RNA can be recovered. Shotgun RNA sequencing of this material allows for simultaneous recovery of viral genomes for transmission tracking as well as understanding of in situ host response [15]. Since the first detection of SARS-CoV-2 in the US in Washington State, the University of Washington Virology Laboratory has performed shotgun RNA sequencing to recover more than 1,500 viral genomes as of August 16, 2020, to understand the evolution and molecular epidemiology of the virus [16,17]. Here, we examine host-specific gene expression differences by SARS-CoV-2 infection status, host age, sex, and viral load in NP swabs from 430 SARS-CoV-2-infected individuals and 54 negative controls.

Results

Between early March and August 16, 2020, the University of Washington Virology Laboratory tested more than 550,000 samples, primarily NP swabs, for infection with SARS-CoV-2 [18]. Thousands of SARS-CoV-2 positive samples, as well as negative controls, have subsequently been metagenomically sequenced, contributing to a detailed understanding of the phylogeny and molecular epidemiology of the virus [16,19]. In this study, we selected a subset of sequenced samples that had sufficient reads pseudoaligned to the human transcriptome (>500,000, median 1.61 × 106) to examine gene expression changes as a result of reverse transcription polymerase chain reaction (RT-PCR)-confirmed SARS-CoV-2 infection. Patient demographics are summarized in Table 1.

[Figure omitted. See PDF.]

Table 1. Patient demographics of SARS-CoV-2 positive and negative samples.

https://doi.org/10.1371/journal.pbio.3000849.t001

SARS-CoV-2 infection induces an antiviral, interferon-mediated host response

We first characterized the genes most differentially expressed (DE) in the nasopharynx as a result of SARS-CoV-2 infection (n = 430 positive, 54 negative). After correcting for batch effects, we found 83 DE genes (padj < 0.1 and absolute log2 fold change > 1) between SARS-CoV-2 positive and negative samples, comprising 41 up-regulated genes and 42 down-regulated genes (S1 Table). Clustering of samples by the 50 most significant DE genes reveals multiple gene expression clusters among SARS-CoV-2 positive samples, while most negative samples cluster together (Fig 1A). Consistent with results from Butler and colleagues [20], SARS-CoV-2 infection induces an interferon-driven antiviral response in the nasopharynx, up-regulating transcripts encoding viral sensors (IFIT1, DDX58), chemokines that attract effector T cells and natural killer (NK) cells (CXCL9, 10, 11), and direct inhibitors of viral replication and function (MX2, RSAD2, HERC5), highlighted in Fig 1A and 1B.

[Figure omitted. See PDF.]

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.

https://doi.org/10.1371/journal.pbio.3000849.g001

To interrogate the global regulatory and signaling programs induced by SARS-CoV-2 infection, we employed Gene Set Enrichment Analysis (GSEA) [21,22] of the 50 Hallmark Gene Sets of the Molecular Signatures Database [23]. Sets with a significant (false discovery rate [FDR] < 0.05) positive enrichment score included Interferon Alpha, Interferon Gamma, and Inflammatory Responses (Fig 1C, S1A Fig). Interestingly, we also found several metabolic pathways negatively enriched, including both Oxidative Phosphorylation and Glycolysis, suggesting a global reduction in production of proteins related to cellular energy production (Fig 1C, S1B Fig). Broad down-regulation of transcripts encoding metabolic machinery may represent either an antiviral response or viral-mediated disruption of host transcripts. We also performed a statistical enrichment test against the Biological Processes Gene Ontology (GO) [24,25]. The most enriched processes (S1C–S1E Fig) are related to either immune responses or translation. In addition to up-regulation of innate antiviral transcripts (Fig 1B, S1D Fig), we also found a consistent down-regulation of transcripts encoding ribosomal proteins (S1E Fig).

ISGs are differentially expressed as a function of viral load

The SARS-CoV-2 receptor ACE2 is an interferon-regulated gene and is up-regulated in response to SARS-CoV-2 infection [8]. We examined the relationship between viral load—defined by the cycle threshold (Ct) of the SARS-CoV-2 nucleocapsid gene region 1 (N1) target during diagnostic PCR—and ACE2. We found that ACE2 expression was associated with increased viral load: median counts of negative, low viral load (N1 Ct > 24), medium viral load (N1 Ct 24–19), and high viral load (N1 Ct < 19) were 0, 1.93, 3.45, and 7.82, respectively (p = 7.46 × 10−13, by Kruskal-Wallis one-way ANOVA; Fig 2A). A similar trend was found for other interferon-induced genes, a subset of which is shown in Fig 2A, including those significantly up-regulated in SARS-CoV-2 infection (CXCL9, OASL, and MX1), negative regulators of inflammation (CD274/PD-L1 and USP18), and monocyte chemoattractant protein-1 (CCL2) [26]. Conversely, the protease required for viral entry, TMPRSS2, was reduced upon viral infection but was not modulated by viral dose, nor were ribosomal proteins (RPL4, RPS6).

[Figure omitted. See PDF.]

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.

https://doi.org/10.1371/journal.pbio.3000849.g002

We next specifically examined gene expression differences between high (N1 Ct < 19, n = 108) relative to low (N1 Ct > 24, n = 99) viral load samples. Fig 2B highlights the 15 most up-regulated and 15 most down-regulated of the 363 total DE genes (adjusted p < 0.1, S2 Table). While genes up-regulated in high viral load samples were dominated by proinflammatory and/or interferon-induced factors such as CXCL9/10, IDO1, and CD80, genes with higher expression in low viral load samples included chemokines for neutrophils (CXCL8, S100A9), and B cell–specific transcripts (FCRL2, IGHG1, IGHM, IGLL5, IGHG2, CD22). Because this suggested differences in immune infiltration as a result of viral load, we performed in silico cell sorting of immune cells using CIBERSORTx [27] and found a higher proportion of transcripts associated with naïve B and T cells, neutrophils, and M2-polarized macrophages in low viral load samples (3.5-, 2.2-, 1.6-, and 1.8-fold increased, respectively), while high viral load samples contained a larger proportion of M1 macrophages, activated NK cells, and activated dendritic cells (DCs) (2.5-, 1.6-, and 1.6-fold up-regulated, respectively; Fig 2C). Levels of transcripts encoding B cell proteins and neutrophil chemokines varied by viral load (Fig 2D).

Human airway epithelial cells recapitulate the cell-intrinsic antiviral response seen in NP swabs

Detection of differential infiltration of antigen-presenting cells and lymphocytes to the nasopharynx in high versus low viral load samples highlights the role that immune cells play in the host response to SARS-CoV-2. To understand whether in vivo infection could be adequately modeled in vitro, we examined gene expression differences in male human airway epithelial (HAE) cells 3 and 7 days post infection with SARS-CoV-2 and compared the DE genes at day 7 to those from SARS-CoV-2 positive versus negative (Fig 1) and high versus low viral load SARS-CoV-2-positive samples (Fig 2). Although HAE are tracheal/bronchial epithelial cells and thus reside lower in the respiratory tract than the nasopharynx, we nevertheless found a consensus set of 19 up-regulated genes that define cell-intrinsic host antiviral responses to SARS-CoV-2 infection (Fig 3A; p = 4.81 × 10−69). When this consensus set was tested for statistical enrichment in the DisGeNET [28] of disease ontologies, we found a high degree of overlap with influenza signature genes (Fig 3B), including a number of interferon-induced genes that mediate the acute antiviral response in the respiratory tract (Fig 3C). Notably, in the HAE cells, there was no sign of induction of an interferon response at 3 days post infection in spite of a 10-fold higher infectious dose of virus used and virus making up 0.3% of reads. At 7 days post infection, SARS-CoV-2 constituted 5.3% of reads.

[Figure omitted. See PDF.]

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−69) 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.

https://doi.org/10.1371/journal.pbio.3000849.g003

Longitudinal analysis of patient-matched samples shows reduced viral load and increased wound healing and humoral immune responses over time

Infection time course may account for the observed differences in immune-related genes in high versus low viral load samples: Patients receiving repeat SARS-CoV-2 testing have a reduced viral load over time [13,29,30]. Although we do not have the ability to tie our data from our large set of positive samples back to the onset of symptoms, we have seen gene expression changes in a small set (n = 3) of matched longitudinal samples with a mean elapsed time of 6.3 days between collections and a mean increase in Ct value of 5.29, representing a 39-fold reduction in viral load (Fig 4A). Enriched GO terms include those related to translation and immune regulation (Fig 4B). Notably, in the second sample collected, we saw increases in genes such as C1QA, -B, and -C and HLA-DQB1 that drive humoral immune responses and those involved in wound healing (APOE, CD36, RHOC), as well as reductions in negative regulators of each process (i.e., TREM1, TFPI). We also tested the 19 SARS-CoV-2 signature genes (Fig 3A) and found reductions in most at the second collection timepoint, although only RSAD2, IFIT2, and HERC5 decreases were statistically significant with 3 samples. We also saw a recovery in the expression of ribosomal proteins over time (Fig 4C). Analysis of data from additional patients with inclusion of more extensive clinical data is currently under investigation in our group.

[Figure omitted. See PDF.]

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: padj < 0.1, white bars: padj > 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.

https://doi.org/10.1371/journal.pbio.3000849.g004

Factors affecting T-cell trafficking are reduced in older adults infected with SARS-CoV-2

Clinically, COVID-19 cases tend to be more severe for older adults and males [1]. No significant difference in N1 Ct was observed based on age or sex (Fig 5A). To understand the differences in host response to SARS-CoV-2 infection, we tested the interaction between infection and age (greater than 60), controlling for non–infection-related age differences with SARS-CoV-2-negative samples. We found only 2 genes altered as a result of the interaction between age and SARS-CoV-2 infection: a 30-fold reduction in production of CXCL11 (Fig 5B), an interferon-induced chemokine for NK and CD8+ T cells, and a 17-fold reduction in polycomb group factor 6 (PCGF6) (S2A Fig), a polycomb repressor complex protein known to play a role in repression of dendritic cell activation [31]. Although we did not find additional genes altered specifically as a result of age in SARS-CoV-2 infection, we did find that CXCL9 and CXCL10 are not induced as strongly in SARS-CoV-2-positive patients aged 60 or older. We also found reduced expression of the receptor for CXCL9/10/11 (CXCR3), the apoptosis-inducing factor GZMB secreted by NK and T cells, and the effector T-cell marker (CD8A). These data suggest that age-related T-cell and NK cell dysfunction [32,33] may play a role in SARS-CoV-2 pathogenesis in older individuals.

[Figure omitted. See PDF.]

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 (padj < 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.

https://doi.org/10.1371/journal.pbio.3000849.g005

Males and females have differences in immune responses to SARS-CoV-2

We performed a similar analysis to evaluate sex differences in SARS-CoV-2 infection and found a total of 19 genes for which the differences in expression based on sex could be attributed to SARS-CoV-2 infection. S2B Fig highlights the top 10 nonredundant enriched GO categories, most of which are related to immune function. In men, we found down-regulation of B cell–specific transcripts (IGHG1 and MSA4A/CD20), down-regulation of the NK cell–activating receptor SLAMF6, and an up-regulation of several inhibitors of NFκB signaling (NDRG1, ARRB2, CD300A, and NFKBID) (Fig 5D). The down-regulation of B cell–specific markers suggests differences in lymphocyte composition and/or trafficking in males. Furthermore, the reduction in NK cell–activating receptors and up-regulation of negative regulators of immune effector function—and resultant throttling of effector function—is consistent with a more severe manifestation of COVID-19 in males.

Discussion

One of the hallmarks of COVID-19 is a dysregulated antiviral immune response. Studies of SARS-CoV, which also employs ACE2/TMPRSS2-mediated entry, have demonstrated that infection does not always result in production of interferons in macrophages and dendritic cell [34], and significantly delayed expression of type II or III interferon in lung cells [10]. Moreover, infection of BALB/c mice with SARS-CoV did not result in detectable interferon beta (IFNβ) until 24 hours, at which point viral titers had nearly reached a peak; lung damage resulting from the subsequent massive infiltration of inflammatory macrophages could be abrogated by pre-treatment with type I interferons [11]. Similar viral kinetics have been observed in SARS-CoV-2-infected patients [13,14] and ferrets [12]. Collectively, these results support a common mechanism by which SARS-CoV and SARS-CoV-2 suppress intracellular viral detection and subsequent interferon induction long enough for viral replication to occur.

Our transcriptomic analysis of NP swabs, which were collected primarily after the onset of symptoms, reveals the robust induction of an interferon response by SARS-CoV-2 infection, similar to that observed by Butler and colleagues [20] and consistent with genes regulated by other coronaviruses [35,36]. The highest levels of individual interferon-responsive genes were seen in samples with the highest viral load. Furthermore, although the in silico cell sorting algorithm we employed has not been specifically validated in NP samples, it detected an enrichment of transcripts associated with inflammatory macrophages and activated DCs and NK cells, 3 primary sources of type I and II interferons, in high viral load samples. When repeat swabs were taken from patients with an average 6.3-day time period between sampling, the interferon response had waned, as had viral load, while transcripts associated with an anti-inflammatory, wound-healing response had increased. Notably, in contrast to the robust expression of type I interferon-regulated transcripts in HAE 7 days after infection with SARS-CoV-2, there was limited evidence of induction of an antiviral response after only 3 days, consistent with a SARS-CoV-like functional repression of intrinsic antiviral immunity and interferon signaling. Clearly, high temporal-resolution transcriptional data from both in vitro culture systems and patient-matched longitudinal samples, as well as clinical data for patients, are necessary to understand the early events in viral pathogenesis and immune-mediated resolution of infection, respectively.

COVID-19 patients frequently develop interleukin-6 (IL-6)-driven cytokine release syndrome (CRS), and elevated serum IL-6 correlates with respiratory failure and poor clinical outcomes [37]. Treatment with the IL-6 receptor blocking antibody tocilizumab has effectively treated COVID-19 symptoms in some patients [38]. We did not see a significant difference in expression of IL-6, nor of other CRS-associated factors such as TNF or VEGF, when we analyzed SARS-CoV-2 positive samples relative to negative, nor in high versus low viral load SARS-CoV-2-positive samples. This could be attributed to the nasopharynx not being a particularly sensitive anatomic location to probe markers of systemic inflammation compared to serum or lower respiratory sites; additional studies will be necessary to ascertain how SARS-CoV-2-induced DE of transcripts in the nasopharynx relates to inflammation and pathology in affected organ systems such the lungs, intestinal tract, and circulatory system. Notably, the dissection of interferon responses in the lower respiratory tract will be important in light of recent findings that type III interferons in particular can disrupt healing and barrier function of lung epithelial cells following viral infection [39,40], a possible counterindication to the proposed use of type III interferon as treatment for COVID-19. Finally, our choice to use a large number of samples at relatively low sequencing depth likely reduced our sensitivity to detect differences in low-abundance and short-lived transcripts like cytokines, including interferons.

One of the more striking patterns we observed is the marked down-regulation of transcription of ribosomal proteins upon SARS-CoV-2 infection (Fig 1B) and the recovery of expression during disease progression (Fig 4C). Global inhibition of host transcription is a strategy employed by many viruses via diverse mechanisms such as disrupting transcriptional pre-initiation complex assembly [41,42] or cleavage of TATA-binding protein [43]. Middle East respiratory syndrome-coronavirus (MERS-CoV) and SARS-CoV nsp1 both cause decay of host mRNA [44,45]; in MERS-CoV, host mRNA degradation results from an endonucleolytic function of nsp1 itself [46]. Nsp1 from SARS-CoV and SARS-CoV-2 share 84% amino acid identity, therefore it is likely that SARS-CoV-2 nsp1 can also function directly or indirectly to promote host RNA degradation. Global down-regulation of host transcription may also be driven in part by SARS-CoV-2 ORF6 protein (which binds to the mRNA export factor RAE1 and nuclear pore protein Nup98 in a similar manner as the vesicular stomatitis virus [VSV] M protein) [47] or the ORF7a protein (which binds to proteins involved in ribosomal assembly and nuclear export) [47].

Finally, understanding age- and sex-related differences in responses to SARS-CoV-2 infection is of critical importance as approximately 90% of SARS-CoV-2 deaths in Washington State have been seen in individuals over age 60. Our data show that in individuals over 60, expression of interferon-induced chemokines is reduced, possibly contributing to a reduction in transcripts for cytotoxic T and NK cells. Immune dysfunctions in older individuals are well-characterized [2,3,32,33] and likely contribute to poorer COVID-19 outcomes; results from clinical trials of type I and III interferons in severely ill patients are likely to further define the role of interferon signaling in older adults [48–50].

Differences in immune responses in males and females are due to a variety of factors, including the effects of sex hormones and the X-linked nature of many immune genes [51]. The bias toward expression of B cell transcripts in females in our study is consistent with higher levels of B cells in females regardless of age [52]. Females also tend to have increased inflammation in response to viral infections [4]. The observed increased expression of inhibitors of nuclear factor kappa-B (NF-κB) in males with SARS-CoV-2 may represent either inappropriate throttling of the antiviral immune response or an adaptive mechanism to reduce deleterious inflammation, a hallmark of COVID-19 pathogenesis.

Collectively, we demonstrate induction of an antiviral response characterized by type I and II interferon induction, which wanes with time and is correlated with viral load. We also find evidence of transcriptional repression by SARS-CoV-2. Lastly, we show that differences in immune responses may underlie disparities in outcomes for 2 higher-risk groups, males and the elderly.

Methods

Ethics statement

Sequencing of excess clinical samples was approved by the University of Washington Institutional Review Board (STUDY00000408) with a consent waiver.

Sample collection, RNA extraction, and RT-PCR

NP swabs of patients with suspected SARS-CoV-2 infection were collected in 3 mL viral transport medium (VTM). Total RNA was extracted from 200 or 140 μL of VTM using either the Roche MagNAPure (Roche, Basel, Switzerland) or Qiagen BioRobot (Qiagen, Hilden, Germany) automated platforms, respectively [53]. RT-PCR for the SARS-CoV2 N1 target was performed on the Applied Biosystems 7500 real time PCR instrument [54,55].

Library preparation and sequencing

Metagenomic next-generation sequencing (mNGS) was performed as previously described [17,56]. Briefly, 18 μL of extracted RNA was treated with Turbo DNAse (ThermoFisher, Waltham, MA). First strand cDNA synthesis was completed using SuperScript IV (ThermoFisher) and random hexamers (Invitrogen, Carlsbad, CA) followed by second strand synthesis by Sequenase version 2.0 (ThermoFisher). The resulting cDNA was purified using either the DNA Clean & Concentrator kit (Zymo, Irvine, California) or 1.6× volumes of AMPure XP beads (Beckman Coulter, Brea, CA). Library preparation was performed using the Nextera XT Kit (Illumina, San Diego, CA). Libraries were cleaned with 0.7× or 0.75× volumes of Ampure beads (Beckman Coulter), quantified using either the Qubit dsDNA HS assay (ThermoFisher) or Quant-iT dsDNA HS assay (ThermoFisher), quality checked by Bioanalyzer or TapeStation (Agilent, Santa Clara, CA), pooled, and sequenced on 1 × 75 bp runs on an Illumina NextSeq or 1 × 101 bp runs on an Illumina NovaSeq.

Pseudoalignment

Raw FASTQ files were adapter and quality trimmed by Trimmomatic version 0.39 [57] using the call “leading 3 trailing 3 slidingwindow:4:15 minlen 20.” Trimmed reads were pseudoaligned to the Ensembl version 96 human transcriptome using Kallisto version 0.46 [58] assuming an average library size of 300 ± 100 bp. Only samples with more than 500,000 pseudoaligned reads were used for RNA sequencing analysis.

DE

Pseudoaligned reads were pre-filtered to remove any genes with average expression of less than one read per sample, then normalized, and DE was calculated with the R package DEseq2 version 1.28.1 [59]. Correction for batch effects was incorporated into the design formula and modeling performed using the Wald test with outlier replacement. Results were deemed significant at a Benjamini-Hochberg adjusted p < 0.1. Gene expression differences attributable to sex or age were incorporated into the design formula as interaction terms.

GSEA

GSEA was performed on normalized counts on GSEA software version 4.0.3 [21,22]. Gene ranking was generated with the Signal2Noise metric and analyzed against the mSigDB Hallmarks version 7.1 gene sets [23].

HAE cultures

The EpiAirway AIR-100 system (MatTek Corporation, Ashland, MA) consists of normal human-derived tracheo/bronchial epithelial cells that have been cultured to form a pseudostratified, highly differentiated mucociliary epithelium closely resembling that of epithelial tissue in vivo. Upon receipt from the manufacturer, HAE cultures were transferred to 6-well plates containing 1.0 mL EpiAirway medium per well (basolateral feeding, with the apical surface remaining exposed to air) and acclimated at 37°C in 5% CO2 for 24 hours prior to experimentation.

Viral growth in HAE

HAE cultures were infected by applying 200 μL of EpiAirway phosphate-buffered saline (MatTeK TEER Buffer) containing 2,000 PFU or 20,000 PFU of infectious clone-derived SARS-CoV-2 expressing a stable mNeonGreen reporter gene (icSARS-CoV-2-mNG) [60] to the apical surface for 90 min at 37°C. At 90 min, the medium containing the inoculum was removed, the apical surface was washed with 200 μL of TEER buffer, and cultures were placed at 37°C. Cultures were fed every other day with 1.0 mL medium via the basolateral side. Media was removed, and cultures were lysed with TRIzol Reagent (ThermoFisher) at 3 days post infection (20,000 PFU challenge) and at 7 days post infection (2,000 PFU challenge). Bam files of viral sequence are deposited in the sequence read archive, NCBI Bioproject PRJNA634194.

HAE RNA sequencing and analysis

RNA from uninfected and infected HAE was extracted using Direct-zol RNA MicroPrep (Zymo). Libraries were generated using the TruSeq Stranded mRNA kit (Illumina) and 2 × 100 bp paired-end reads sequenced on a Novaseq. Pseudoalignment using Kallisto version 0.44 and DE analysis was performed as earlier.

Statistics and visualization

All calculations were performed in R version 4.0.0. Statistical enrichment tests of GO [24,25] and DisGeNET [28] pathways were performed in the clusterProfiler R package [61]. Exact test of the intersection of HAE and NP swab consensus genes was performed with the R package SuperExactTest [62]. Images were generated using packages including DOSE [63], ggplot2, pheatmap, and VennDiagram.

Supporting information

S1 Fig. DE gene sets and GO Biological Process terms in SARS-CoV-2 NP swabs.

(A) Enrichment plots of gene sets significantly (FDR < 0.05) positively enriched in SARS-CoV-2 samples. (B) Enrichment plots of gene sets significantly (FDR < 0.05) negatively enriched in SARS-CoV-2 samples. (C) Top 20 Biological Process GO terms for which DE genes in SARS-CoV-2 samples are overrepresented. “Number Enriched” is the number of SARS-CoV-2 DE genes that belong to each GO term. (D) Fold change of genes belonging to GO term “defense response to virus.” (E) Fold change of genes belonging to GO term “SRP-dependent cotranslational protein targeting to membrane.” Raw data available in the GEO Repository, accession GSE152075.

https://doi.org/10.1371/journal.pbio.3000849.s001

(TIF)

S2 Fig. Age and sex differences in gene expression upon SARS-CoV-2 expression.

(A) 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. (B) Top 10 Biological Process GO terms in which genes defining the male versus female response to virus are overrepresented. Raw data available in the GEO Repository, accession GSE152075.

https://doi.org/10.1371/journal.pbio.3000849.s002

(TIF)

S1 Table. DE genes in SARS-CoV-2-positive samples relative to negative.

https://doi.org/10.1371/journal.pbio.3000849.s003

(CSV)

S2 Table. DE genes in SARS-CoV-2 high viral load samples relative to low viral load.

https://doi.org/10.1371/journal.pbio.3000849.s004

(CSV)

Acknowledgments

The authors thank Amin Addetia and Joshua Lieberman for their thoughtful comments to improve the manuscript.

Citation: Lieberman NAP, Peddu V, Xie H, Shrestha L, Huang M-L, Mears MC, et al. (2020) In vivo antiviral host transcriptional response to SARS-CoV-2 by viral load, sex, and age. PLoS Biol 18(9): e3000849. https://doi.org/10.1371/journal.pbio.3000849

References

1. OpenSAFELY: factors associated with COVID-19-related hospital death in the linked electronic health records of 17 million adult NHS patients. | medRxiv [Internet]. [cited 2020 Jun 5]. Available from: https://www.medrxiv.org/content/10.1101/2020.05.06.20092999v1.

2. Mueller AL, McNamara MS, Sinclair DA. Why does COVID-19 disproportionately affect older people? Aging. 2020 May 29;12(10):9959–81. pmid:32470948

3. McElhaney JE, Verschoor CP, Andrew MK, Haynes L, Kuchel GA, Pawelec G. The immune response to influenza in older humans: beyond immune senescence. Immun Ageing. 2020 May 7;17(1):10.

4. vom Steeg LG, Klein SL. SeXX Matters in Infectious Disease Pathogenesis. PLoS Pathog. 2016 Feb;12(2):e1005374. pmid:26891052

5. Channappanavar R, Fett C, Mack M, Ten Eyck PP, Meyerholz DK, Perlman S. Sex-based differences in susceptibility to SARS-CoV infection. J Immunol Baltim Md 1950. 2017 May 15;198(10):4046–53.

6. Meng Y, Wu P, Lu W, Liu K, Ma K, Huang L, et al. Sex-specific clinical characteristics and prognosis of coronavirus disease-19 infection in Wuhan, China: A retrospective study of 168 severe patients. PLoS Pathog. 2020 Apr 28;16(4):e1008520. pmid:32343745

7. Hoffmann M, Kleine-Weber H, Schroeder S, Krüger N, Herrler T, Erichsen S, et al. SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor. Cell. 2020 Apr 16;181(2):271–280.e8. pmid:32142651

8. Ziegler CGK, Allon SJ, Nyquist SK, Mbano IM, Miao VN, Tzouanas CN, et al. SARS-CoV-2 Receptor ACE2 Is an Interferon-Stimulated Gene in Human Airway Epithelial Cells and Is Detected in Specific Cell Subsets across Tissues. Cell. 2020 May;181(5):1016–1035.e19. pmid:32413319

9. Hou YJ, Okuda K, Edwards CE, Martinez DR, Asakura T, Dinnon KH, et al. SARS-CoV-2 Reverse Genetics Reveals a Variable Infection Gradient in the Respiratory Tract. Cell. 2020 May;S0092867420306759.

10. Menachery VD, Eisfeld AJ, Schäfer A, Josset L, Sims AC, Proll S, et al. Pathogenic influenza viruses and coronaviruses utilize similar and contrasting approaches to control interferon-stimulated gene responses. mBio. 2014 May 20;5(3):e01174–01114. pmid:24846384

11. Channappanavar R, Fehr AR, Vijay R, Mack M, Zhao J, Meyerholz DK, et al. Dysregulated Type I Interferon and Inflammatory Monocyte-Macrophage Responses Cause Lethal Pneumonia in SARS-CoV-Infected Mice. Cell Host Microbe. 2016 Feb 10;19(2):181–93. pmid:26867177

12. Blanco-Melo D, Nilsson-Payant BE, Liu W-C, Uhl S, Hoagland D, Møller R, et al. Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19. Cell. 2020 28;181(5):1036–1045.e9. pmid:32416070

13. He X, Lau EHY, Wu P, Deng X, Wang J, Hao X, et al. Temporal dynamics in viral shedding and transmissibility of COVID-19. Nat Med. 2020 May;26(5):672–5. pmid:32296168

14. Zou L, Ruan F, Huang M, Liang L, Huang H, Hong Z, et al. SARS-CoV-2 Viral Load in Upper Respiratory Specimens of Infected Patients. N Engl J Med. 2020 Mar 19;382(12):1177–9. pmid:32074444

15. Greninger AL, Chen EC, Sittler T, Scheinerman A, Roubinian N, Yu G, et al. A metagenomic analysis of pandemic influenza A (2009 H1N1) infection in patients from North America. PLoS ONE. 2010 Oct 18;5(10):e13381. pmid:20976137

16. Fauver JR, Petrone ME, Hodcroft EB, Shioda K, Ehrlich HY, Watts AG, et al. Coast-to-Coast Spread of SARS-CoV-2 during the Early Epidemic in the United States. Cell. 2020 28;181(5):990–996.e5. pmid:32386545

17. Peddu V, Shean RC, Xie H, Shrestha L, Perchetti GA, Minot SS, et al. Metagenomic analysis reveals clinical SARS-CoV-2 infection and bacterial or viral superinfection and colonization. Clin Chem. 2020 May 7.

18. Randhawa AK, Fisher LH, Greninger AL, Li SS, Andriesen J, Corey L, et al. Changes in SARS-CoV-2 Positivity Rate in Outpatients in Seattle and Washington State, March 1-April 16, 2020. JAMA. 2020 May 8.

19. Bedford T, Greninger AL, Roychoudhury P, Starita LM, Famulare M, Huang M-L, et al. Cryptic transmission of SARS-CoV-2 in Washington State. medRxiv. 2020 Apr 16;2020.04.02.20051417.

20. Butler DJ, Mozsary C, Meydan C, Danko D, Foox J, Rosiene J, et al. Shotgun Transcriptome and Isothermal Profiling of SARS-CoV-2 Infection Reveals Unique Host Responses, Viral Diversification, and Drug Interactions. bioRxiv. 2020 Apr 20;2020.04.20.048066.

21. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005 Oct 25;102(43):15545–50. pmid:16199517

22. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003 Jul;34(3):267–73. pmid:12808457

23. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015 Dec 23;1(6):417–25. pmid:26771021

24. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000 May;25(1):25–9. pmid:10802651

25. The Gene Ontology Consortium. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res. 2019 08;47(D1):D330–8. pmid:30395331

26. Hubel P, Urban C, Bergant V, Schneider WM, Knauer B, Stukalov A, et al. A protein-interaction network of interferon-stimulated genes extends the innate immune system landscape. Nat Immunol. 2019 Apr;20(4):493–502. pmid:30833792

27. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019 Jul;37(7):773–82. pmid:31061481

28. Piñero J, Ramírez-Anguita JM, Saüch-Pitarch J, Ronzano F, Centeno E, Sanz F, et al. The DisGeNET knowledge platform for disease genomics: 2019 update. Nucleic Acids Res. 2020 08;48(D1):D845–55. pmid:31680165

29. Wölfel R, Corman VM, Guggemos W, Seilmaier M, Zange S, Müller MA, et al. Virological assessment of hospitalized patients with COVID-2019. Nature. 2020;581(7809):465–9. pmid:32235945

30. Liu Y, Yan L-M, Wan L, Xiang T-X, Le A, Liu J-M, et al. Viral dynamics in mild and severe cases of COVID-19. Lancet Infect Dis. 2020;20(6):656–7. pmid:32199493

31. Boukhaled GM, Cordeiro B, Deblois G, Dimitrov V, Bailey SD, Holowka T, et al. The Transcriptional Repressor Polycomb Group Factor 6, PCGF6, Negatively Regulates Dendritic Cell Activation and Promotes Quiescence. Cell Rep. 2016 16;16(7):1829–37. pmid:27498878

32. Hazeldine J, Lord JM. The impact of ageing on natural killer cell function and potential consequences for health in older adults. Ageing Res Rev. 2013 Sep;12(4):1069–78. pmid:23660515

33. Li M, Yao D, Zeng X, Kasakovski D, Zhang Y, Chen S, et al. Age related human T cell subset evolution and senescence. Immun Ageing. 2019 Sep 11;16(1):24.

34. Frieman M, Baric R. Mechanisms of Severe Acute Respiratory Syndrome Pathogenesis and Innate Immunomodulation. Microbiol Mol Biol Rev MMBR. 2008 Dec;72(4):672–85. pmid:19052324

35. Fung TS, Liu DX. Human Coronavirus: Host-Pathogen Interaction. Annu Rev Microbiol. 2019 08;73:529–57. pmid:31226023

36. Kindler E, Thiel V, Weber F. Interaction of SARS and MERS Coronaviruses with the Antiviral Interferon Response. Adv Virus Res. 2016;96:219–43. pmid:27712625

37. Chen G, Wu D, Guo W, Cao Y, Huang D, Wang H, et al. Clinical and immunological features of severe and moderate coronavirus disease 2019. J Clin Invest. 2020 May 1;130(5):2620–9. pmid:32217835

38. Xu X, Han M, Li T, Sun W, Wang D, Fu B, et al. Effective treatment of severe COVID-19 patients with tocilizumab. Proc Natl Acad Sci. 2020 May 19;117(20):10970–5. pmid:32350134

39. Major J, Crotta S, Llorian M, McCabe TM, Gad HH, Priestnall SL, et al. Type I and III interferons disrupt lung epithelial repair during recovery from viral infection. Science. 2020 Jun 11.

40. Broggi A, Ghosh S, Sposito B, Spreafico R, Balzarini F, Lo Cascio A, et al. Type III interferons disrupt the lung epithelial barrier upon viral recognition. Science. 2020 Jun 11.

41. Dasgupta A, Scovell WM. TFIIA abrogates the effects of inhibition by HMGB1 but not E1A during the early stages of assembly of the transcriptional preinitiation complex. Biochim Biophys Acta. 2003 Jun 19;1627(2–3):101–10. pmid:12818428

42. Di Valentin E, Bontems S, Habran L, Jolois O, Markine-Goriaynoff N, Vanderplasschen A, et al. Varicella-zoster virus IE63 protein represses the basal transcription machinery by disorganizing the pre-initiation complex. Biol Chem. 2005 Mar;386(3):255–67. pmid:15843171

43. Kundu P, Raychaudhuri S, Tsai W, Dasgupta A. Shutoff of RNA polymerase II transcription by poliovirus involves 3C protease-mediated cleavage of the TATA-binding protein at an alternative site: incomplete shutoff of transcription interferes with efficient viral replication. J Virol. 2005 Aug;79(15):9702–13. pmid:16014932

44. Lokugamage KG, Narayanan K, Nakagawa K, Terasaki K, Ramirez SI, Tseng C-TK, et al. Middle East Respiratory Syndrome Coronavirus nsp1 Inhibits Host Gene Expression by Selectively Targeting mRNAs Transcribed in the Nucleus while Sparing mRNAs of Cytoplasmic Origin. J Virol. 2015 Nov 1;89(21):10970–81. pmid:26311885

45. Kamitani W, Narayanan K, Huang C, Lokugamage K, Ikegami T, Ito N, et al. Severe acute respiratory syndrome coronavirus nsp1 protein suppresses host gene expression by promoting host mRNA degradation. Proc Natl Acad Sci. 2006 Aug 22;103(34):12885–90. pmid:16912115

46. Nakagawa K, Narayanan K, Wada M, Popov VL, Cajimat M, Baric RS, et al. The Endonucleolytic RNA Cleavage Function of nsp1 of Middle East Respiratory Syndrome Coronavirus Promotes the Production of Infectious Virus Particles in Specific Human Cell Lines. J Virol. 2018 01;92(21).

47. Gordon DE, Jang GM, Bouhaddou M, Xu J, Obernier K, White KM, et al. A SARS-CoV-2 protein interaction map reveals targets for drug repurposing. Nature. 2020 Apr 30;1–13.

48. Park A, Iwasaki A. Type I and Type III Interferons–Induction, Signaling, Evasion, and Application to Combat COVID-19. Cell Host Microbe [Internet]. 2020 May 27 [cited 2020 Jun 7]; Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7255347/.

49. Fragkou PC, Belhadi D, Peiffer-Smadja N, Moschopoulos CD, Lescure F-X, Janocha H, et al. Review of trials currently testing treatment and prevention of COVID-19. Clin Microbiol Infect [Internet]. 2020 May 23 [cited 2020 Jun 7]; Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7245266/.

50. Hung IF-N, Lung K-C, Tso EY-K, Liu R, Chung TW-H, Chu M-Y, et al. Triple combination of interferon beta-1b, lopinavir–ritonavir, and ribavirin in the treatment of patients admitted to hospital with COVID-19: an open-label, randomised, phase 2 trial. The Lancet. 2020 May 30;395(10238):1695–704.

51. Klein SL, Flanagan KL. Sex differences in immune responses. Nat Rev Immunol. 2016;16(10):626–38. pmid:27546235

52. Teixeira D, Longo-Maugeri IM, Santos JLF, Duarte YAO, Lebrão ML, Bueno V. Evaluation of lymphocyte levels in a random sample of 218 elderly individuals from São Paulo city. Rev Bras Hematol E Hemoter. 2011;33(5):367–71.

53. Lieberman JA, Pepper G, Naccache SN, Huang M-L, Jerome KR, Greninger AL. Comparison of Commercially Available and Laboratory Developed Assays for in vitro Detection of SARS-CoV-2 in Clinical Laboratories. J Clin Microbiol. 2020 Apr 29;

54. Nalla AK, Casto AM, Huang M-LW, Perchetti GA, Sampoleo R, Shrestha L, et al. Comparative Performance of SARS-CoV-2 Detection Assays Using Seven Different Primer-Probe Sets and One Assay Kit. J Clin Microbiol. 2020 26;58(6).

55. Perchetti GA, Nalla AK, Huang M-L, Zhu H, Wei Y, Stensland L, et al. Validation of SARS-CoV-2 detection across multiple specimen types. J Clin Virol Off Publ Pan Am Soc Clin Virol. 2020 13;128:104438.

56. Greninger AL, Zerr DM, Qin X, Adler AL, Sampoleo R, Kuypers JM, et al. Rapid Metagenomic Next-Generation Sequencing during an Investigation of Hospital-Acquired Human Parainfluenza Virus 3 Infections. J Clin Microbiol. 2017;55(1):177–82. pmid:27795347

57. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinforma Oxf Engl. 2014 Aug 1;30(15):2114–20.

58. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7. pmid:27043002

59. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. pmid:25516281

60. Xie X, Muruato A, Lokugamage KG, Narayanan K, Zhang X, Zou J, et al. An Infectious cDNA Clone of SARS-CoV-2. Cell Host Microbe. 2020 13;27(5):841–848.e3. pmid:32289263

61. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics J Integr Biol. 2012 May;16(5):284–7.

62. Wang M, Zhao Y, Zhang B. Efficient Test and Visualization of Multi-Set Intersections. Sci Rep. 2015 Nov 25;5:16923. pmid:26603754

63. Yu G, Wang L-G, Yan G-R, He Q-Y. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinforma Oxf Engl. 2015 Feb 15;31(4):608–9.