Initial replication of SARS-CoV-2 in the upper respiratory tract is required to establish infection, and the replication level correlates with the likelihood of viral transmission. Here, we examined the role of host innate immune defenses in restricting early SARS-CoV-2 infection using transcriptomics and biomarker-based tracking in serial patient nasopharyngeal samples and experiments with airway epithelial organoids. SARS-CoV-2 initially replicated exponentially, with a doubling time of ∼6 h, and induced interferon-stimulated genes (ISGs) in the upper respiratory tract, which rose with viral replication and peaked just as viral load began to decline. Rhinovirus infection before SARS-CoV-2 exposure accelerated ISG responses and prevented SARS-CoV-2 replication. Conversely, blocking ISG induction during SARS-CoV-2 infection enhanced viral replication from a low infectious dose. These results show that the activity of ISG-mediated defenses at the time of SARS-CoV-2 exposure impacts infection progression and that the heterologous antiviral response induced by a different virus can protect against SARS-CoV-2.
SARS-CoV-2 enters the body and first replicates in the upper respiratory tract, usually peaking within the first week of infection. The level of viral replication in the nasopharynx correlates strongly with the likelihood of transmission, with most transmission occurring in the first week (Cevik et al., 2021; He et al., 2020). Furthermore, significant viral replication following SARS-CoV-2 exposure is likely a prerequisite (although certainly not the only factor) for COVID-19 progression. Thus, successful host responses at this early stage could reduce viral transmission and progression to COVID-19 disease.
Innate immune responses such as the antiviral IFN response are often critical for curtailing early viral replication at mucosal surfaces; however, recent work shows that SARS-CoV-2 antagonizes and delays this response (Banerjee et al., 2020b; Blanco-Melo et al., 2020; Konno et al., 2020; Xia et al., 2020; Zhou et al., 2020; Martin-Sancho et al., 2021; Ravindra et al., 2021). Nonetheless, there is strong evidence that the IFN response protects against severe COVID-19, from both genetic and autoantibody studies, and from in vitro studies showing protection of cells using recombinant IFNs. However, it is unclear how early in infection these defenses function (Lee and Shin, 2020; Meffre and Iwasaki, 2020; Zhang et al., 2020; Bastard et al., 2020; Vanderheiden et al., 2020; Pairo-Castineira et al., 2021; Lokugamage et al., 2020). A recent study linked rapid induction of virus-specific T cells to early viral clearance (Tan et al., 2021). However, the timing and role of the IFN response during initial SARS-CoV-2 replication in the upper respiratory tract has not been reported. This is due in part to the difficulty of obtaining samples from the start of infection, when the infection is often asymptomatic, and also due to the lack of established, practical methods to track these responses in human subjects.
The mucosal IFN response is initiated when pattern recognition receptors within epithelial cells and immune cells sense general features shared by many viruses, such as common structural features of viral RNA. This recognition event triggers expression of type I and type III IFNs and IFN-stimulated genes (ISGs). Secreted IFNs, in turn, bind to cell surface receptors on nearby cells, amplifying ISG expression and creating an antiviral state in the mucosal barrier (Schneider et al., 2014; Odendall and Kagan, 2015; Iwasaki, 2012). Most respiratory viruses induce an IFN response in the airway, although the kinetics and magnitude vary, and most are also highly susceptible to this host defense mechanism if it is sufficiently activated during active viral replication.
Recent work on the host response to SARS-CoV-2 in the airway epithelium or model cell lines has shown that SARS-CoV-2 can trigger IFN and ISG expression via the viral RNA sensor MDA5, albeit with delayed kinetics, but that this response in some cases was ineffective at curtailing viral replication; furthermore, while pretreatment with exogenous IFN blocked SARS-CoV-2 infection, IFN was much less effective if added after the infection was established (Yin et al., 2021; Rebendenne et al., 2021; Hsin et al., 2021,Preprint). Single-cell analysis of epithelial cells from the nasopharynx of COVID-19 patients showed ISG induction in epithelial cells primarily from patients with mild-to-moderate disease, but not severe disease, consistent with the observation that patients with severe disease are more likely to have genetic or acquired deficiencies in IFN signaling pathways (Ziegler et al., 2021 Preprint). Together, these findings indicate that epithelial cells are capable of mounting an IFN response to SARS-CoV-2 but that the timing and magnitude of the IFN response relative to viral replication may be critical for determining the course of infection.
In this study, we sought to capture early host–virus dynamics in the human nasopharynx using serial patient samples and study the functional consequences of modulating the host innate immune response in primary human airway epithelial air–liquid interface (ALI; organoid) cultures. While SARS-CoV-2 is known to antagonize and delay the IFN response, we hypothesized that this defense mechanism might still be critical for restricting initial SARS-CoV-2 replication under some circumstances, depending on the relative kinetics of viral replication and ISG induction.
Due to SARS-CoV-2 screening and testing practices at our hospital at the start of the pandemic in March 2020, we were able to obtain nasopharyngeal (NP) swab samples from patients at multiple time points after infection, including serial samples collected close to the start of infection in asymptomatic subjects. Using transcriptomic and biomarker-based analysis of these samples, we observed robust ISG induction in the airway mucosa in diverse COVID-19 patients. Longitudinal analysis revealed a dynamic mucosal response, with a rise and fall in ISG expression tracking with the rise and fall in viral load. Using the organoid model, we found that prior infection with rhinovirus (RV) accelerated ISG activity early in SARS-CoV-2 infection and completely prevented SARS-CoV-2 replication, an effect that was reversed by preventing ISG induction. Conversely, blocking homologous ISG induction during SARS-CoV-2 infection enhanced the replication rate in the setting of a low infectious dose. Together, these results demonstrate that innate immunity in the airway of human subjects is highly dynamic and demonstrate that innate immune responses profoundly influence initial infection with SARS-CoV-2.
SARS-CoV-2 induces an IFN response in the nasopharynx across diverse patient groups
The host response to SARS-CoV-2 in the upper respiratory tract and its relationship to viral replication at the start of infection are not well defined. To initially characterize host responses during SARS-CoV-2 infection in vivo, we performed RNA sequencing (RNA-seq) on NP swab RNA from SARS-CoV-2–positive patients (n = 30) and SARS-CoV-2–negative healthcare worker controls (n = 8). Patients included outpatients and patients admitted to the hospital of both sexes, who ranged in age from 20 to 90 yr, with the majority being >60 yr of age (Table S1). Samples varied in viral load over five orders of magnitude as assessed by quantitative RT-PCR (RT-qPCR) for the SARS-CoV-2 N1 gene or by read mapping of NP RNA to the SARS-CoV-2 genome, with a strong correlation seen between PCR Ct value for the SARS-CoV-2 N1 gene and viral RNA reads by RNA-seq (r2 = 0.8380, P < 0.0001).
Of reads mapping to the human genome, 1,770 RNAs differed significantly between SARS-CoV-2–positive patients and control subjects (Fig. 1 A). These included 1,567 protein-coding genes, of which 1,245 (79.4%) were enriched and 322 reduced in patients relative to controls. The most significantly enriched genes in the nasopharynx of SARS-CoV-2 patients were known ISGs, including three NP ISG transcripts previously shown by our group to accurately identify patients with viral respiratory infection (OASL, IFIT2, and CXCL10; Fig. 1 A and Table S2; Landry and Foxman, 2018). Analysis of Ingenuity pathways and transcription factor–binding sites associated with enriched transcripts demonstrated activation of multiple pathways related to ISG induction in SARS-CoV-2 patients compared with controls, as well as other pathways linked to innate immunity, leukocyte recruitment, and initiation of mucosal inflammatory responses (Fig. 1, B–D; and Table S3). Notably, transcripts regulated by neuronal tissue–associated transcription factors were enriched in control subjects compared with SARS-CoV-2–positive patients (Fig. S1). Many of these transcripts were previously shown to be more highly expressed in the human olfactory epithelium and/or olfactory sensory neurons compared with the respiratory epithelium, consistent with reports that SARS-CoV-2 targets the olfactory epithelial cells and possibly neurons in the nasopharynx (Olender et al., 2016; Fodoulian et al., 2020).
Examination of gene expression across patient samples revealed several patterns (Fig. 1, E–G). First, the 45 most significantly enriched genes were all ISGs, according to the Interferome database (Rusinova et al., 2013) and consistent with previous reports (Mick et al., 2020; Lieberman et al., 2020). ISGs appeared to be coregulated within individual patients (i.e., patients with high expression of one ISG tended to have high expression of other ISGs; Fig. 1 E). This was also demonstrated by analysis of the correlation between reads for different ISGs across samples (Fig. 1 F). Second, ISG expression appeared to be loosely correlated with viral load, with those patients with the highest viral load (Fig. 1 E, left) tending to have higher ISG expression than those with the lowest viral loads (Fig. 1 E, right). However, while all SARS-CoV-2–positive samples showed enrichment of ISG expression compared with controls, there was no clear relationship between ISG expression level and sex, age, or outpatient/admitted status in this diverse set of 30 samples. Notably, IFNλ1 and IFNγ, but not other IFN subtypes, showed significant enrichment in COVID-19 patients compared with controls, although expression levels were relatively low compared with many ISGs, similar to findings in a single-cell RNA-seq (scRNA-seq) study of epithelial cells from the nasopharynx of COVID-19 patients (Table S2; Ziegler et al., 2021, Preprint). Direct comparison of DEGs in patient groups with distinct clinical characteristics (male versus female, older versus younger, and inpatient versus outpatient) showed no significant differences in ISG expression, although there were some other notable biological differences in NP gene expression that have not been reported previously. Male and female patients showed differential NP expression of a subset of sex-linked genes, including the innate immune associated RNA helicase DDX3Y (Table S4; Szappanos et al., 2018). Also, younger patients (<55 yr) showed enrichment for mRNAs expressed by the repair-associated Hillock cells of the airway epithelium (KRT4 and KRT13) compared with older patients (>70 yr), and a similar pattern was seen in outpatients compared with admitted patients (Table S5 and Table S6) In addition, some previously reported patterns were also noted, such as increased expression of a subset of airway mucins in COVID-19 patients compared with controls (Table S2; Biering et al., 2021 Preprint).
NP CXCL10 protein correlates with ISG RNA in NP swabs
Next, we sought a quantitative biomarker of ISG induction in order to evaluate this response across a larger set of samples and over time in individual patients. We previously showed that NP CXCL10 is detected in the viral transport medium during other acute viral respiratory infections and correlates with expression of ISGs at the mRNA level (Landry and Foxman, 2018). Therefore, we measured the level of CXCL10 protein in the NP swab–associated viral transport medium using ELISA. We observed a significant positive correlation between NP CXCL10 protein level with the NP mRNA level of CXCL10 in SARS-CoV-2–positive patient NP samples (Fig. 1, E and G). Together, these results indicated that across subjects with diverse clinical presentations, SARS-CoV-2 induced a robust IFN response in the nasopharynx and that the NP CXCL10 protein level correlated with ISG expression at the RNA level.
NP host response correlates with SARS-CoV-2 viral load across patient groups
Using CXCL10 as a biomarker, we evaluated a larger set of samples from SARS-CoV-2–positive patients tested in our healthcare system in March and April of 2020 (n = 140) in order to gain further insight into the relationships among viral replication, disease status, and host response to infection (Fig. 2; and Fig. S2, A and B). First, we examined NP CXCL10 level in SARS-CoV-2–positive individuals who were tested as outpatients and not admitted to the hospital compared with those who were admitted. Notably, we observed significantly higher NP CXCL10 levels in outpatients compared with admitted patients (Fig. 2 A; P = 0.0019). To understand the reason for this, we first examined patient age, since as a group, admitted patients were significantly older than outpatients (16 yr older on average; Fig. S2 C). However, there was no correlation between age and CXCL10 level (Fig. S2 D). Next, we examined viral load. We were initially surprised to find that admitted patients had significantly lower viral loads than outpatients (Fig. 2 B). This suggested that the main factor driving CXCL10 level was viral load. Supporting this idea, correlation analysis showed a significant positive correlation between NP CXCL10 level and viral load by RT-qPCR (for all patients, r2 = 0.2030, P < 0.0001; Fig. 2 C). This correlation was also seen in separate analyses of outpatients and admitted patients, but there was no significant difference in the slope of the CXCL10 vs. viral load correlation between these groups, although there was a trend toward a higher slope in outpatients (Fig. 2 C). We also observed no significant relationship between sex and NP CXCL10 or sex and viral load in this sample set (Fig. S2, E–H).
Prior work on SARS-CoV-2 has shown that the NP viral load is highest in the first few days of infection and that the more severe symptoms of COVID-19 requiring hospitalization occur in the second or third week of infection (Cevik et al., 2021). Therefore, we hypothesized that admitted patients may have had lower viral loads at the time of testing than outpatients because they presented later in infection, after peak viral replication in the nasopharynx. Consistently, outpatients tended to report fewer days of symptoms before testing compared with admitted patients, although this information was only available for a subset of patients (approximately one third, n = 44; Fig. 2 D).
NP CXCL10 reveals kinetics of host response to SARS-CoV-2 at the start of infection
To further evaluate the relationship between viral replication and the innate antiviral response in the nasopharynx, we examined viral load and NP CXCL10 data in longitudinal patient samples. First, we examined viral load from 29 inpatients from March 12 and April 30, 2020 for whom we had at least eight sequential test results for SARS-CoV-2 with at least the first sample tested using the CDC assay (our clinical laboratory also had other testing platforms; Fig. 3 and Table S7). In March and April 2020, most serial testing in our healthcare system was aimed at patient clearance for discharge. Consistently, the majority of patients (15/29) showed low viral loads (Ct N1 > 21), which remained low throughout the time course. Another common pattern was high viral load in the first sample (Ct N1 < 20) followed by a decline in viral load over time, similar to patterns reported in the literature for patients who presented close to the start of symptomatic illness (7/29 patients; Fig. 3 A). These patients showed high CXCL10 level in the sample with peak viral load and a decline in NP CXCL10 after the viral load had decreased (Fig. 3 B). One patient had consistently high NP viral loads for 20 d and did not survive (not shown). Finally, a third pattern was seen in several patients (6/29), in which the first sample had a low viral load that subsequently increased to a high peak level (Ct N1 <20) and then decreased over time (Fig. 3, C–H). This pattern is consistent with patient presentation close to the start of infection. Two of these patients had no symptoms of SARS-CoV-2 at presentation, and the virus was detected incidentally on screening during hospitalization for other reasons (Fig. 3, C and D; and Table S7). One of these patients had an inconclusive test result and a positive test result on the same day, 12 h later (first test, N1 not detected, N2 Ct = 38.4; second test, N1 Ct = 34.6, N2 Ct = 35.4), suggesting that this might have been the first day of infection for this patient (L2, Fig. 3 C). The other four patients presented with acute symptoms, including fever and, in some cases, cough and/or shortness of breath (Fig. 3, E–H; and Table S7). NP CXCL10 in these patients was undetectable or low in the first positive sample with low viral load, rose with viral RNA, and then subsequently declined as viral load declined. Together, the longitudinal data show a correlation with NP CXCL10 and viral load in individual patients over time, similar to what we observed across 140 patients tested at a single time point (Fig. 2 C). Notably, for patient L2, the only patient for whom three samples were available before peak viral load, there appeared to be a delay between CXCL10 production relative to viral replication during the first few days of infection, suggesting that initially, viral replication outpaced the host innate immune response in the nasopharynx.
SARS-CoV-2 RNA increases exponentially, with a doubling time of ∼6 h, and infection induces a robust ISG response in organoids by day 4 after infection
To further evaluate the kinetics of SARS-CoV-2 replication and host response, we performed a time course of SARS-CoV-2 infection using primary human airway epithelial cells grown at ALI, which differentiate into organoids with beating cilia and mucus production, recapitulating the airway mucosal surface in vivo. Reminiscent of what we observed during SARS-CoV-2 infection in vivo (Fig. 3, C–H), intracellular viral load in the organoids increased rapidly between 1 and 72 h and then plateaued between 72 and 96 h; in the apical wash, viral load increased exponentially starting slightly later, from 24 to 96 h after infection (Fig. 4, A and B). ISGs were also induced, with ISG mRNA levels and CXCL10 protein level in the basolateral medium increasing markedly from 72 to 96 h (Fig. 4, C–F). Notably, a very high level of CXCL10 protein was produced by infected epithelia (∼4 ng/ml by 96 h), consistent with the strong NP CXCL10 signal observed in the nasopharynx in vivo (Fig. 4 F). While in theory the CXCL10 measured in frozen NP swabs could reflect intracellular CXCL10, this result indicates that epithelial cells likely secrete CXCL10, consistent with our prior study in which CXCL10 was detected in the NP swab–associated media from virus-positive patients when cells were first removed from fresh samples by centrifugation (Landry and Foxman, 2018). Next, we checked for IFNλ1 at the mRNA and protein level, since the RNA-seq data showed significant induction of IFNλ1 and IFNγ in the nasopharynx of COVID-19 patients and IFNλ1 is known to be highly produced by epithelial cells, whereas IFNγ is usually produced by immune cells. We observed robust induction of IFNλ1 mRNA and a small but significant increase in IFNλ1 protein by 96 h after infection (Fig. 4, G and H).
Viral replication in organoid cultures appeared to follow exponential growth during the first few days of infection. Although viruses do not replicate by doubling per se, an exponential increase in viral load indicates that during this phase of growth, the increase in virus at a given time is proportional to the current amount. In this situation, the doubling time, which describes the rate of exponential increase, is a useful mathematical descriptor for viral replication rate. Therefore, we used curve fitting to exponential growth to estimate the doubling time of SARS-CoV-2 in organoids, which was 5.858 h (95% confidence interval [CI], 4.85–7.357 h; Fig. 4 I) and in the apical wash, which was 9 h (95% CI, 7.438–11.45 h; Fig. 4 J).
Although we did not evaluate viral load by plaque assay, other studies that measured SARS-CoV-2 titer during replication in similar culture systems also observed a period of exponential replication at the start of infection, consistent with our results. However, not unexpectedly, the reported kinetics vary with experimental conditions, such as viral load in the inoculum and incubation temperature (Hou et al., 2020; Zhu et al., 2020; V’kovski et al., 2021). Thus, we sought to estimate the replication kinetics in vivo using serial patient samples. For patient L2, viral load from the first three SARS-CoV-2–positive time points also appeared to follow exponential growth; therefore, we used the same method to estimate the SARS-CoV-2 doubling time in vivo from these data, which was 6.454 h (95% CI, 4.261–13.30 h based on 3 y-values; Fig. 4 K). For all other patients for whom viral load increased in serial samples (Fig. 3), we had only one sample before peak viral load. We asked what the doubling times for SARS-CoV-2 would be in these samples if we assumed exponential replication between the first and peak viral RNA values. The calculated doubling times across patients ranged from 3.048 to 6.509 h for samples ≤5 d apart (Fig. 4 L). For the two patients with a larger sampling interval (L44, L12), calculated doubling times were 9.455 and 12.58 h, although these calculations would be expected to overestimate doubling time if the second sample was taken after viral replication had plateaued or began to decline (Fig. 4 L). Together, the in vitro and in vivo results indicate that SARS-CoV-2 replicates exponentially during the first few days after infection before the peak host antiviral response, with an average doubling time of ∼6 h.
RV infection increases mRNA level of dACE2, but not full-length ACE2, in an IRF3-dependent manner
In organoid cultures, SARS-CoV-2 appeared to induce ISGs with delayed timing relative to viral replication, consistent with other cellular infection models (Yin et al., 2021). We considered that other physiological exposures, such as a prior infection with a different virus, could potentially accelerate these ISG responses. We focused on RV, since it is the most frequently detected virus in the human upper respiratory tract (Turner, 2007). Based on previous studies, we expected that RV might curtail infection by enhancing antiviral responses but also potentially promote infection by increasing expression of SARS-CoV-2 entry receptors (Wu et al., 2020; Ziegler et al., 2020). Using conditions previously established to produce robust viral amplification, we infected organoids with RV (HRV-01A, multiplicity of infection [MOI] 0.05) and evaluated the effect on expression of ACE2, the SARS-CoV-2 entry receptor. ACE2 was originally reported to be an ISG, but a subsequent study reported that full-length ACE2, which functions as an entry receptor, is not an ISG, and that a truncated form, dACE2, is an ISG but is not a functional SARS-CoV-2 entry receptor (Onabajo et al., 2020; Ziegler et al., 2020). Consistent with this finding, we observed that dACE2 was significantly induced by RV infection (∼10-fold). Induction was prevented by BX795, a signaling inhibitor that prevents IRF3 nuclear translocation and that we previously showed prevents ISG induction by RV in airway epithelial cells (Clark et al., 2009; Wu et al., 2020; Fig. S3). Notably, full-length ACE2 expression was also significantly increased by RV infection, albeit to a much lesser extent than dACE2 (approximately twofold), and this change was not abrogated by BX795 (Fig. S3) RV infection had no effect on expression of TMPRSS2, a cofactor for SARS-CoV-2 entry (not shown). Together, these results indicate that at the mRNA level, RV induces mRNA for dACE2 ∼10-fold by an IRF3-dependent mechanism, as expected for an ISG (Onabajo et al., 2020), and induces a small increase in mRNA for the SARS-CoV-2 entry receptor, full-length ACE2, approximately twofold by a different mechanism.
Prior RV infection blocks SARS-CoV-2 replication
Next, we evaluated SARS-CoV-2 replication and ISG induction following infection of airway epithelial organoids, with or without prior RV infection (Fig. 5 A). RV induced expression of ISGs by day 3 after infection (Fig. 5 B). Importantly, RV-infected organoids appeared to be healthy at 3 d after infection, with cultures showing beating cilia similar to uninfected cultures (Video 1 and Video 2). SARS-CoV-2 viral load increased exponentially in infected cultures without prior RV infection, as observed previously (Fig. 4), but showed essentially no increase between 24 and 72 h when cultures had been exposed to RV 3 d prior (Fig. 5 C). Nonlinear regression analysis was consistent with exponential replication of SARS-CoV-2 from 24–72 h after infection in mock-pretreated cultures, a significant difference from RV-preinfected cultures (Fig. 5 C, P value). These results indicate that despite the slight increase in mRNA for full-length ACE2 triggered by RV (Fig. S3), the net result of RV infection 3 d before SARS-CoV-2 infection is to block replication of SARS-CoV-2. Next, we performed IFNλ1 and IFNβ ELISA on the basolateral media from organoid cultures infected with SARS-CoV-2 only or RV followed by SARS-CoV-2. For IFNλ1, we observed minimal IFNλ1 production in cultures infected with SARS-CoV-2 only but high levels of IFNλ1 from organoids sequentially infected with RV followed by SARS-CoV-2 (Fig. 5 D). There was no significant increase in IFNβ in SARS-CoV-2–infected or RV/SARS-CoV-2–coinfected cultures (Fig. 5 E), consistent with the literature showing that type III IFNs, not type I, are the major IFNs produced by respiratory epithelial cells in response to viral infection (Khaitov et al., 2009; Okabayashi et al., 2011). Consistently, evaluation of ISG expression over the course of infection showed that at early time points of SARS-CoV-2 infection, ISGs were significantly more highly expressed in RV-preinfected cultures that in cultures infected with SARS-CoV-2 without prior RV exposure (Fig. 5 F). This included several ISGs that have been previously reported to limit coronavirus replication or for which polymorphisms are linked to disease severity of SARS-CoV or SARS-CoV-2, including ISG15, BST2 (tetherin), LY6E, and OAS1-3 (Hamano et al., 2005; He et al., 2006; Ma et al., 2014; Taylor et al., 2015; Pfaender et al., 2020; Martin-Sancho et al., 2021; Pairo-Castineira et al., 2021; Fig. 5 F).
RV induces an IFN response in bystander cells
To better understand the timing and breadth of the epithelial host response to RV that appeared to limit SARS-CoV-2 replication, we evaluated ISG expression over time for 5 d and examined ISG expression and viral infection at the single-cell level. Time-course analysis showed that following inoculation, RV replicated robustly, peaking at 24 h after infection, and then declined significantly but was still detectable by RT-qPCR at day 5, a time point corresponding to 48 h after SARS-CoV-2 infection in the sequential infection experiment (Fig. S4 A). ISG expression increased and decreased in parallel with viral replication but was still significantly higher than in mock-treated cells at day 5 after RV infection (Fig. S4, B–D). In addition, a high level of IFNλ1 protein was detected in the basolateral media of RV-infected cultures at day 5 after infection (Fig. S4 E). Next, we performed scRNA-seq of mock-treated and RV-infected organoid cultures to evaluate the host ISG response at single-cell resolution at day 5 after infection At this time point, RV viral RNA reads were detected in only 70 out of 4,200 cells (Fig. 5 G). The infected cells were predominantly ciliated cells but included all major cell types, consistent with the HRV-01A entry receptor low-density lipoprotein receptor being ubiquitously highly expressed throughout the culture (Fig. S4, F and G). Although RV RNA was only detected in a small subset of cells in infected cultures at day 5 (1.67%), ISGs were elevated in all cells compared with mock-treated cultures (Fig. 5 G), indicating that RV infection induces a robust antiviral response in uninfected cells and consistent with detection of a high level of IFNλ1 in the basolateral media (Fig. S4 E). Although it is possible that some of the ISG-expressing cells in the day 5 culture had been previously been infected with RV, RV is generally considered to be a lytic infection and to infect only a small proportion of cells within the airway mucosa; thus, it is likely that the decline in viral load from day 1 to day 5 reflects death of infected cells and the majority of cells in the culture at day 5 are bystander cells (Jacobs et al., 2013). Consistent with this idea, RV infection did not significantly alter the cellular composition of organoid cultures on day 5 (Fig. S4 H). Together, these lines of evidence indicate that RV infection induces robust ISG expression in all cells of infected cultures, many of which are bystander cells.
RV blocks SARS-CoV-2 replication through heterologous ISG induction
Next, to test whether suppression of SARS-CoV-2 replication in the presence of RV was dependent upon the host innate immune response, we blocked ISG induction by including the signaling inhibitor BX795 throughout single or sequential infections (Fig. 6 A and Fig. S5). Notably, under these conditions (MOI 0.5, 72 h), the homologous ISG response to SARS-CoV-2 did not appear to be limiting viral replication, since BX795 suppressed ISG induction by SARS-CoV-2 but did not alter SARS-CoV-2 viral load within the culture or apical wash (Figs. 6 C and S5). However, BX795 dramatically altered SARS-CoV-2 replication in the context of RV infection. As seen previously (Fig. 5), prior infection with RV suppressed SARS-CoV-2 replication by >1,000-fold (Fig. 6 C). However, following BX795 pretreatment, SARS-CoV-2 replication was restored (Fig. 6 C). These results indicate that ISG induction is required for the inhibition of SARS-CoV-2 replication by prior RV infection.
RV RNA was detected at much lower levels than SARS-CoV-2 RNA at this time point (72 h after SARS-CoV-2 infection, 6 d after RV infection), and showed a slight reduction during SARS-CoV-2 coinfection without BX795, but significantly higher levels when both viruses were present in the presence of BX795 (Fig. 6 D). This result indicates that the antiviral response also limits RV replication in the setting of SARS-CoV-2 coinfection at this time point. In other words, in the presence of an intact antiviral response, the viral load of both viruses was reduced, but if ISG induction was inhibited, the viral load was equal (for SARS-CoV-2) or higher (for RV) during coinfection. This experiment models sequential viral infections in a host with an intact IFN response (both viruses decrease during coinfection) compared with a host with a deficient IFN response (equal or greater replication of both viruses during coinfection).
Homologous ISG response limits SARS-CoV-2 replication in a low-MOI infection
Next, to probe whether there are any conditions under which the homologous host response to SARS-CoV-2 might limit early viral replication, we tested the effect of BX795 using a 10-fold-lower level of virus in the inoculum (MOI 0.05), since under these conditions, SARS-CoV-2 replication would potentially be more sensitive to suppression by ISGs (Fig. 6, E–H). In this experiment, BX795 treatment led to an ∼10-fold increase in organoid-associated virus and an ∼300-fold increase in virus shedding into the apical wash at 72 h after infection. A trend toward an increase (5–10 times) was also seen at 96 h after infection, although due to variability among replicates, this difference was not statistically significant (Fig. 6, F and G). Next, based on the increase in viral RNA from 1 to 72 h during the low-MOI infection, we estimated the effect of BX795 on SARS-CoV-2 doubling time in organoid culture. Assuming exponential growth between 1 and 72 h, the SARS-CoV-2 doubling time was 5.127 h without BX795 (95% CI, 3.889–7.518 h) and 3.578 h with BX795 (95% CI, 3.499–3.661 h; Fig. 6 H). These results indicate that homologous ISG induction during the first 3 d of infection limits SARS-CoV-2 replication rate in the setting of a low infectious dose.
Taken together, the data from patient samples and organoid culture experiments show that SARS-CoV-2 undergoes exponential replication at the start of infection and induces ISGs that peak just before the fall in NP viral load. The heterologous response to a previous RV infection accelerated ISG induction and inhibited SARS-CoV-2 replication in its target tissue, the airway epithelium. Furthermore, the homologous ISG response induced by the virus itself significantly curtailed viral replication from a low starting viral load. These results demonstrate that ISG-mediated defenses in the airway epithelium are dynamic and define specific instances in which this host defense profoundly limits SARS-CoV-2 replication at the start of infection.
SARS-CoV-2 transmission fuels the continuation of the COVID-19 pandemic. Peak viral replication in the upper respiratory tract occurs during the first few days of infection and correlates directly with transmission (Cevik et al., 2021; He et al., 2020). However, host–virus interactions at this stage of disease are not fully understood, in part due to the difficulties in capturing samples from this early stage of the infection and in part due to lack of straightforward methods to measure such responses. Here, we establish such methods and show that SARS-CoV-2 replicates exponentially at the start of infection and robustly induces ISGs in the nasopharynx within the first week, but not from the start, of infection. Through follow-on mechanistic studies in an organoid model, we show that varying the timing of ISG induction profoundly impacts the extent of viral replication (Figs. 5, 6, and S5).
Due to the nature of exponential growth, even a small change in the replication rate, as we observed with inhibition of ISG induction during low-MOI infection (Fig. 6, E–H), could have a profound impact on peak viral load. For example, if conditions allowed the viral doubling time to decrease by 2 h from 6 h to 4 h, this would lead to a 64-fold greater NP viral load after 72 h. NP viral load has been shown to correlate with viral transmission, an issue that has come into focus recently due to the emergence of SARS-CoV-2 variants with enhanced transmission (Cevik et al., 2021; He et al., 2020; Singanayagam et al., 2020; Centers for Disease Control and Prevention, 2021a). Our results suggest that viral mutations that enhance IFN/ISG antagonism would enable a faster doubling time and compel studying at this aspect of the biology as a possible mechanism for increased transmissibility in emerging viral variants.
Conversely, factors that enhance IFN-mediated defenses would be expected to reduce peak viral load, viral transmission, and host susceptibility. In this study, we focused on heterologous innate immunity induced by RV, the most frequent cause of common colds (Turner, 2007). Community studies have revealed that RV is much more prevalent in the upper respiratory tract than previously appreciated, detected in ∼15% of asymptomatic subjects at a given time and often at much higher rates in young children (Foxman and Iwasaki, 2011; Jartti et al., 2008; Byington et al., 2015). RVs and other respiratory viruses have also been shown to induce ISGs in the nasopharynx in vivo, in both the presence and absence of symptoms (Landry and Foxman, 2018; Wolsk et al., 2016; Yahya et al., 2017; Yu et al., 2019). Therefore, heterologous ISG induction may be important in defense against closely spaced viral infections, as we observed in a previous study of RV and influenza A (Wu et al., 2020). However, the outcome of closely spaced infections could also be influenced by virus-specific features, including ISG antagonism and entry receptor usage. The predicted effect of prior RV infection on SARS-CoV-2 replication was unclear, since RV induces transcription of the SARS-CoV-2 entry receptor, full-length ACE2, in the airway epithelium (Fig. S3) Also, SARS-CoV-2 has been shown to antagonize both ISG induction and activity at multiple levels (Banerjee et al., 2020a; Sa Ribero et al., 2020; Xia et al., 2020; Martin-Sancho et al., 2021). However, our experiments showed that prior RV infection protected against replication of SARS-CoV-2, that this protection was dependent upon ISG induction, and that a significant bystander IFN response was detected in cells throughout the epithelium for at least 5 d after RV infection, even after the RV itself was largely cleared (Figs. 5, 6, and S4). These results indicate that the protective effect of heterologous ISG induction throughout the epithelium by RV predominated over other effects, potently suppressing SARS-CoV-2 replication.
While host–pathogen interactions are often studied one at a time in experimental models, in the human upper respiratory tract, exposure to multiple microbes simultaneously or in series is a frequent occurrence and likely an important influence on innate immunity. The idea that viral interference could be shaping broad patterns of respiratory virus transmission has only recently begun to gain traction, in part due to epidemiological data suggesting interference among RNA respiratory viruses (Isaacs and Lindenmann, 1957; Greer et al., 2009; Linde et al., 2009; Casalegno et al., 2010; Ånestad and Nordbø, 2011; Schultz-Cherry, 2015; Karppinen et al., 2016; Nickbakhsh et al., 2019; Wu et al., 2020). Heterologous ISG induction could be particularly important for host defense against a virus like SARS-CoV-2, since it would preempt mechanisms SARS-CoV-2 has in place to antagonize and delay IFN and ISG induction in response to its own replication.
The extent to which heterologous ISG responses have affected and will affect SARS-CoV-2 transmission during the pandemic remains unclear. The SARS-CoV-2 pandemic and the public health response have disrupted the status quo in many ways, including dramatically reducing the circulation of common respiratory viruses (Olsen et al., 2020; Jones, 2020; Cowling et al., 2020; Yeoh et al., 2020; Sullivan et al., 2020; Centers for Disease Control and Prevention, 2021b). This change has led to speculation that diminished population adaptive immunity will lead to a surge of infections when these viruses again recirculate, and thus far, early data from school reopenings indicate a likely robust resurgence of RV (Baker et al., 2020; Fong et al., 2021; Poole et al., 2020). Our findings suggest that heterologous innate immune responses could be a mitigating factor, since the model predicts an upper limit in the extent to which IFN-sensitive viruses can simultaneously or sequentially infect the same host (see model, Fig. 7 A). This would also potentially slow viral transmission in the population by reducing the number of susceptible hosts at any given time (Fig. 7 B). Importantly, we are proposing this model based on our data to highlight the importance of considering viral interference in epidemiological patterns that may emerge in the coming year, but we are not arguing that potential effects of heterologous innate immunity outweigh the proven benefits of public health measures for directly preventing SARS-CoV-2 transmission.
There are several important caveats to our study. First, the organoid model contains only epithelial cells, but homologous IFN responses induced by SARS-CoV-2 during initial viral replication are likely amplified by immune cells during natural infection. Based on the early and robust induction of CXCL10, we would expect leukocytes that express the cognate receptor CXCR3, such as T cells and natural killer cells, to be recruited and participate in antiviral defense; thus, in vivo innate responses could be significantly amplified compared with what we observe in vitro. Second, attenuation of the second infection is not the only possible outcome of a viral coinfection. We studied RV, which in most individuals results in asymptomatic infection or common colds. However, for viruses that damage the lung, such as influenza, residual tissue damage could exacerbate rather than protect against a subsequent viral respiratory illness. Studying virus–host–virus interactions for different virus pairs in physiological models is an important future direction of this study. Finally, the concept of protection based on heterologous induction of IFN responses assumes a host with intact innate immune defenses, but this is not always the case. Both genetic and environmental variables can prevent robust ISG induction in response to RV and other respiratory viruses (Lamborn et al., 2017; Asgari et al., 2017; Zhang et al., 2020; Foxman et al., 2015; Mihaylova et al., 2018; Kudo et al., 2019). In our experimental model, there was a profound difference in the outcome of RV-SARS-CoV-2 sequential infection in the presence and absence of an intact host-cell IFN response. With an intact host response, viral loads of both viruses decreased, but when the IFN response was blocked, viral loads of both viruses were equal to or higher than in single infections (Figs. 6 and S5) This result illustrates that the expected outcome of sequential viral infections is likely to be very dependent upon host innate immune status. Other factors include the extent to which coinfecting viruses trigger or block ISG induction and the kinetics of ISG induction relative to the interval between infections.
In sum, our results demonstrate an important role for IFN-mediated defenses in curtailing SARS-CoV-2 replication at the start of infection, including heterologous innate immune responses induced by prior RV infection. These results, and our findings in longitudinal patient samples, support the concept that airway innate immunity is dynamic, with innate immune defense rapidly changing in response to current and recent viral infections. Our findings also demonstrate that ISG-mediated defenses can profoundly curtail SARS-CoV-2 replication under certain conditions and compel further studies of the role of heterologous innate immunity in protecting against SARS-CoV-2 and other respiratory viruses.
Materials and methods
Experimental model and subject details
Human experimental guidelines approval statement
The use of clinical samples and data in this study was approved by the Yale University’s Human Investigation Committee (institutional review board protocol ID #200002765). Procedures for testing residual clinical samples and recording linked patient data, followed by sample and data deidentification, were evaluated, and the requirement for specific patient consent was waived. In vitro experiments used primary human cells obtained from Lonza Bioscience. Lonza Bioscience guarantees that all tissue used for human cell products is ethically obtained with donor informed consent in accordance with processes approved by an institutional review board or comparable independent review body.
We used residual NP samples remaining after clinical testing for CXCL10 measurements and transcriptome analysis. Swab-associated viral transport medium was stored at −80°C following clinical testing and thawed just before ELISA assay or RNA isolation for RNA-seq. Clinical information, including age, sex, virology results, and specific features of clinical course (including presenting symptoms, hospital admission, and length of stay), was extracted from the electronic medical record and recorded, after which samples were assigned a study code and deidentified. In the clinical laboratory, SARS-CoV-2 was detected in most samples using an Emergency Use Authorization–approved TaqMan assay detecting the CDC targets N1, N2, and RNaseP (Centers for Disease Control and Prevention, 2020). In some longitudinal samples, SARS-CoV-2 was diagnosed with the commercial Cepheid assay (Cepheid, 2021); in this case, RT-qPCR for the CDC N1 gene was repeated using RT-qPCR TaqMan assay for the CDC N1 gene as described previously (catalog no. 10006600; Integrated DNA Technologies; Vogels et al., 2020).
Primary human bronchial epithelial cells
Primary human bronchial epithelial cells from healthy adult donors were obtained commercially (Lonza Bioscience) and cultured at ALI according to the manufacturer’s instructions using reduced hydrocortisone (Stem Cell Technologies). Cells were allowed to differentiate for 4 wk, by which time they displayed beating cilia and mucus production. This project used cells from eight different healthy adult donors (four male and four female). Each result was confirmed with independent experiments using organoids grown from at least two different cell donors.
RV 1A (HRV-01A; VR-481; ATCC) was amplified in H1-HeLa cells (CRL-1985; ATCC) and used to generate a high-titer filtered cell lysate. Viral titer was determined by plaque assay as reported previously (Foxman et al., 2015). SARS-CoV-2 (USA-WA1/2020; BEI Resources) was generously provided by the Wilen laboratory (Yale University, New Haven, CT). Virus was cultured on Vero E6 cells, a filtered supernatant was used as the virus stock, and titer was determined by plaque assay as described previously (Ravindra et al., 2021). We confirmed that the cell lysate/supernatant did not contain significant levels of IFN or other molecules that activate ISG expression by mock-inoculating BEAS2B cells with UV-inactivated stocks. This did not induce ISG expression after 24-h incubation, in contrast to robust ISG induction in response to recombinant IFNβ.
RNA isolation from clinical samples
At the time of accessioning, the residual viral transport medium from clinical samples was stored at −80°C. Upon thawing, RNA was isolated from 140 µl transport medium using the Qiagen Viral RNA isolation kit per the manufacturer’s instructions (reference 52904; Qiagen), and one aliquot was reserved for ELISA.
Library preparation and RNA-seq
RNA samples were quantified and checked for quality using the Agilent 2100 Bioanalyzer Pico RNA Assay. Library preparation was performed using Kapa Biosystems Kapa HyperPrep Kit with RiboErase (human, rat, mouse) in which samples were normalized with a total RNA input of 25 ng. Libraries were amplified using 15 PCR cycles. Libraries were validated using Agilent TapeStation 4200 D1000 assay and quantified using the Kapa Library Quantification Kit for Illumina Platforms kit. Libraries were diluted to 1.3 nM and pooled at 1.25% each of an Illumina NovaSeq 6000 S4 flowcell using the XP workflow to generate 25 million read pairs/sample.
RNA-seq data analysis
Low-quality reads were trimmed and adaptor contamination was removed using Trim Galore (v0.5.0, https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Trimmed reads were mapped to the human reference genome (hg38) using HISAT2 (v2.1.0; Kim et al., 2019). Gene expression levels were quantified using StringTie (v1.3.3b) with gene models (v27) from the GENCODE project (Pertea et al., 2015). Differentially expressed genes (adjusted P value < 0.05, fold-change [FC] cutoff = 2) were identified using DESeq2 (v1.22.1; Love et al., 2014). To avoid the unexpected outlier replacement for sex-linked genes, we turned off the outlier replacement option in the male versus female comparison by setting minReplicatesForReplace = Inf for the DESeq() function in the DESeq2 package.
Visualization of RNA-seq data
Protein-coding genes differentially expressed in SARS-CoV-2–positive versus negative control were visualized on a volcano plot, with an x axis cutoff l log2FC = 10. All differentially expressed RNAs are included in Table S1 (n = 1,770). Significantly differentially expressed transcripts were defined as those with log2FC > 1 and adjusted P value < 0.05. Heatmap shows gene expression levels of top 45 most significant DEGs using minimum-to-maximum scaling of normalized read counts. Pathway analysis was performed using Ingenuity Pathway Analysis (version 01–16). Transcription factor motif enrichment analysis was performed using Cytoscape (version 3.8.1) with the iRegulon plug-in (version 1.3; Janky et al., 2014).
In vitro infections
We infected primary human bronchial epithelial cells differentiated at ALI with HRV-01A, SARS-CoV-2, or both. For SARS-CoV-2, high-MOI infection was MOI 0.5, and low-MOI infection was MOI 0.05. For HRV-01A, MOI 0.1 was used, as this was the minimum viral inoculum that reproducibly led to robust HRV-01A viral replication in ALI cultures based on prior studies.
To evaluate the effect of RV on subsequent infection with SARS-CoV-2, we infected with each virus individually or sequentially and examined the time course of viral amplification and ISG induction. To formally test whether prior exposure to RV inhibits SARS-CoV-2 replication through activation of the host-cell IFN response, we performed sequential infection studies in the presence of BX795. For BX795 inhibitor experiments, 6 µM BX795 (Sigma-Aldrich) was added to the media 1 d before the first infection. On the day of the first infection, the medium was replaced with fresh media with or without 6 µM BX795. On the day of the second infection, 150 μl additional media was added to the basolateral chamber, with or without 6 µM BX795. Results shown are representative of at least two independent experiments using primary bronchial epithelial cells from different donors, with four to six biological replicates per condition in each experiment.
For RT-qPCR, RNA was isolated from each well of differentiated epithelial cells using the QIAGEN RNeasy kit by incubating each 24-well insert with 350 µl lysis buffer at room temperature for 5 min, followed by RNA isolation and cDNA synthesis using iScript cDNA synthesis kit (Bio-Rad). To quantify viral RNA and mRNA levels for ISGs and the housekeeping gene HPRT, qPCR was performed using SYBR green iTaq universal (Bio-Rad) per the manufacturer’s instructions. Viral RNA was quantified using primers to the RV genome. Viral RNA per ng total RNA is graphed as FC from the limit of detection (40 cycles of PCR) as 240-Ct. ISG mRNA levels are graphed as FC from mock-treated cells or are presented relative to the level of mRNA for the housekeeping gene HRPT (2−ΔΔCT). RT-qPCR for SARS-CoV-2 within cultures was performed using the previously described TaqMan assay for the CDC N1 gene with primers and probes provided by Integrated DNA Technologies (catalog no. 10006600; Integrated DNA Technologies). RT-qPCR for SARS-CoV-2 in apical was performed using a combined reverse transcription and qPCR reaction using the Luna Universal Probe One-Step RT-qPCR Kit (New England Biolabs). RT-qPCR for ISG expression during SARS-CoV-2 infection was assessed with both HPRT and β-actin housekeeping gene controls to confirm that the HPRT level was stable following viral infection.
The following primers were used for RT-qPCR with SYBR green: HPRT (forward: 5′-TGGTCAGGCAGTATAATCCAAAG-3′; reverse: 5′-TTTCAAATCCAACAAAGTCTGGC-3′), β-actin (forward: 5′-CCTGGCACCCAGCACAAT-3′; reverse: 5′-GCCGATCCACACGGAGTACT-3′), ISG15 (forward: 5′-CATCTTTGCCAGTACAGGAGC-3′; reverse: 5′-GGGACACCTGGAATTCGTTG-3′), RSAD2 (forward: 5′-TCGCTATCTCCTGTGACAGC-3′; reverse: 5′-CACCACCTCCTCAGCTTTTG-3′), MX1 (forward: 5′-AGAGAAGGTGAGAAGCTGATCC-3′; reverse: 5′-TTCTTCCAGCTCCTTCTCTCTG-3′), IFITM3 (forward: 5′-ATCGTCATCCCAGTGCTGAT-3′; reverse: 5′-ATGGAAGTTGGAGTACGTGG-3′), IFIT2 (forward: 5′-CCTCAAAGGGCAAAACGAGG-3′; reverse: 5′-CTGATTTCTGCCTGGTCAGC-3′), CXCL10 (forward: 5′-CCTGCAAGCCAATTTTGTCC-3′; reverse: 5′-ATGGCCTTCGATTCTGGATTC-3′), LY6E (forward: 5′-GCATTGGGAATCTCGTGACA-3′; reverse: 5′-ATGGAAGCCACACCAACATT-3′), BST2 (forward: 5′-CACACTGTGATGGCCCTAAT-3′; reverse: 5′-TGTAGTGATCTCTCCCTCAAGC-3′), IFITM3 (forward: 5′-ATCGTCATCCCAGTGCTGAT-3′; reverse: 5′-ATGGAAGTTGGAGTACGTGG-3′), OAS1 (forward: 5′-GCTCCTACCCTGTGTGTGTGT-3′; reverse: 5′-TGGTGAGAGGACTGAGGAAGA-3′), OAS3 (forward: 5′-GAAAACTGTCAAGGGAGGCTC-3′; OAS3 reverse: 5′-CCCTCTGGTCCACATAGCTC-3′), HRV-01A (forward: 5′-CAGGCCAAATTAAAGTCAATAAGC-3′; reverse: 5′-AGGCTGAAGTTTGGTTTTGC-3′), Fl-ACE2 (forward: 5′-GGGCGACTTCAGGATCCTTAT-3′; reverse: 5′-GGATATGCCCCATCTCATGATG-3′), and d-ACE2 (forward: 5′-GGAAGCAGGCTGGGACAAA-3′; reverse: 5′-AGCTGTCAGGAAGTCGTCCATT-3′).
Quantification and statistical analysis and RT-qPCR data
GraphPad Prism (version 9.0.0) was used for data analysis. Data are shown as mean ± SEM. Statistical significance of differences between conditions was determined by t tests (two tailed) for patient data and by Mann–Whitney tests for in vitro experiments. Linear regression analysis was used to determine association between clinical parameters, such as viral load and NP CXCL10 in clinical samples, and test the null hypothesis that the slope of the association was significantly different from zero. Nonlinear regression analysis was used to fit viral growth to an exponential curve (exponential growth with log(population)) to determine virus doubling times and test the null hypothesis that one curve fit both datasets for the SARS-CoV-2 growth curve with and without RV preinfection. P < 0.05 was considered statistically significant.
CXCL10 levels in cell-free residual NP swab samples or basolateral media of organoid cultures were quantified using a solid-phase sandwich ELISA (catalog no. DY266; R&D Systems). Briefly, frozen viral transport medium from residual NP swab samples or basolateral medium from organoid cultures was thawed on ice and centrifuged to remove cell debris. ELISA was performed according to the manufacturer’s instructions. IFNλ1 protein concentration in culture basolateral media was measured using a commercial IFNλ1 ELISA kit (catalog no. DY7246; R&D Systems) or IFNβ ELISA kit (DY814-05; R&D Systems). Basolateral media from mock-infected and infected cultures was collected from the bottom chamber of the ALI culture, UV inactivated, and stored at −80°C. The ELISA assay was performed on neat and diluted media per the manufacturer’s instructions.
scRNA-seq of ALI organoid cultures
Library preparation and sequencing
Organoid cultures were digested with trypsin/EDTA to form a single-cell suspension. The 10X Genomics single-cell 3′ protocol was used to produce Illumina-ready sequencing libraries with standard Illumina paired-end constructs according to the manufacturer’s instructions.
Analysis of scRNA-seq data
All downstream analyses were implemented using R version 3.6.3 and the package Seurat v3.1.4 (Stuart et al., 2019). Gene expression matrix of mock-infected and RV1A-infected samples was first individually analyzed in this procedure. Genes expressed in fewer than three cells and cells expressing fewer than 200 genes were discarded. A distribution histogram of unique molecular identifier (UMI) count in all cells was made, and cells with <8,000 UMI counts were discarded. This resulted in a matrix of 21,086 genes expressing in 4,511 cells in the mock sample and a matrix of 21,195 genes expressing in 4,200 cells in the RV1A-infected sample. The raw counts were normalized using the Seurat function NormalizeData with normalization.method = “LogNormalize” and scale.facto = 10,000. All genes were scaled using Seurat function ScaleData with default parameters. Variable features were determined using method “vst.” The top 2,000 variable features were used for principal-component analysis. Graph-based clustering was performed individually on mock- and RV1A-infected samples. The K-nearest neighbor graph was built by FindNeighbors function with first 20 principal components and k.param = 10. Louvain clustering was done using FindCluster function with resolution = 0.8 for both samples. Clusters were separately annotated in the mock- and RV1A-infected samples using the following markers of the major cell groups in the airway epithelium: basal cycling (MKI67 and HIST1H4C), basal (KRT14, KRT15, and KRT5), Hillock (SPINK5 and KRT13), secretory (SCGB1A1, BPIFA1, and LYPD2), ciliated (CAPS, PIFO, and MORN2), ionocyte (FOXI1), and pulmonary neuroendocrine cells (markers include AZGP1 and AVIL; Plasschaert et al., 2018). Each cluster found by clustering was assigned to one of the above seven major groups. A group of developing ciliated cells with marker cyclin-O was found in the mock sample and was merged into the ciliated cell cluster. After cluster annotation, mock and infected samples were merged together to produce t-distributed stochastic neighbor embedding (tSNE) maps and make comparisons. Both samples used same normalization method and the gene expression level was rescaled. Dimensionality reduction was performed by latent semantic analysis using Seurat function RunLSI with the first 50 singular values. The tSNE maps were then produced with the first 50 dimensions and perplexity 30. The color coding of tSNE plots used cell type, sample source, viral read per cell and expression levels of genes of interest.
Online supplemental material
Fig. S1 shows neuron-associated transcripts that are enriched in the nasopharynx of healthy controls compared with SARS-CoV-2–positive patients. Fig. S2 shows the relationship between NP CXCL10 levels in SARS-CoV-2–positive patients and age, sex, and inpatient/outpatient status, showing a significant association between NP CXCL10 and outpatient versus admitted status but no significant associations between NP CXCL10 and age or sex within these groups. Fig. S3 shows that mRNA for dACE2, a truncated form of ACE2, is induced ∼10-fold by RV in an IRF3-dependent manner and that full-length ACE2, the SARS-CoV-2 entry receptor, is also induced, but to a lesser extent (approximately twofold), and the induction is not IRF3 dependent. Fig. S4 shows the time course of RV infection and associated IFN and ISG induction in organoid cultures over the course of 5 d. Fig. S5 shows the viral load shed from the apical surface of organoid cultures during high-MOI SARS-CoV-2 infection and ISG expression with and without BX795. Table S1 and Table S7 describe clinical characteristics of patients described in Fig. 1 and Fig. 3, respectively. Table S2, Table S3, Table S4, Table S5, and Table S6 provide information about differentially expressed genes in SARS-CoV-2 patients compared with controls, as shown in Fig. 1. Video 1 and Video 2 show cilia movement in organoid cultures at day 3 after mock inoculation or infection with RV.
scRNA-seq data are available in the GEO database under accession no. GSE164982. RNA-seq data derived from patient samples have been deposited to dbGaP under accession no. phs002433.v1.p1. To access these data, users may apply for access to the dbGaP data repository (https://dbgap.ncbi.nlm.nih.gov/aa/wga.cgi?login=&page=login).
We thank Maureen Owen, Greta Edelman, the entire staff of the Yale New Haven Hospital Clinical Virology laboratory, and Amy Likens for their dedicated assistance. We thank Craig Wilen and Wilen laboratory members for valuable help and advice and providing SARS-CoV-2 virus stock. We also thank Bryan Pasqualucci and Christopher Castaldi at the Yale Center for Genomic Analysis.
Funding was provided by Fast Grants for COVID-19 from Emergent Ventures at the Mercatus Center of George Mason University (E.F. Foxman); the Yale Department of Laboratory Medicine and COVID-19 Dean’s Fund (E.F. Foxman); the China Scholarship Council Yale World Scholars Fellowship (B. Wang); the Gruber Foundation Fellowship (T.A. Watkins); and the National Institutes of Health (grant T32AI007019 to T.A. Watkins).
Author contributions: Conceptualization, N.R. Cheemarla and E.F. Foxman; Data Curation, N.R. Cheemarla, T.A. Watkins, and B. Wang; Methodology, N.R. Cheemarla, T.A. Watkins, V.T. Mihaylova, B. Wang, D. Zhao, and E.F. Foxman; Investigation, N.R. Cheemarla, T.A. Watkins, V.T. Mihaylova, and B. Wang; Formal Analysis, B. Wang, D. Zhao, and G. Wang; Resources, M.L. Landry; Supervision, E.F. Foxman; Visualization: N.R. Cheemarla, B. Wang, and E.F. Foxman; Validation, N.R. Cheemarla and E.F. Foxman; Writing – Original Draft, E.F. Foxman; Writing – Review & Editing, N.R. Cheemarla, T.A. Watkins, V.T. Mihaylova, B. Wang, D. Zhao, G. Wang, M.L. Landry, and E.F. Foxman; Funding Acquisition, E.F. Foxman.
Disclosures: E.F. Foxman reported a patent application to WO2019/217296 A1 pending, and E.F. Foxman and M. L. Landry reported a patent application to WO2018/071498 pending. No other disclosures were reported.