Single cell resolution of SARS-CoV-2 tropism, antiviral responses, and susceptibility to therapies in primary human airway epithelium

The human airway epithelium is the initial site of SARS-CoV-2 infection. We used flow cytometry and single cell RNA-sequencing to understand how the heterogeneity of this diverse cell population contributes to elements of viral tropism and pathogenesis, antiviral immunity, and treatment response to remdesivir. We found that, while a variety of epithelial cell types are susceptible to infection, ciliated cells are the predominant cell target of SARS-CoV-2. The host protease TMPRSS2 was required for infection of these cells. Importantly, remdesivir treatment effectively inhibited viral replication across cell types, and blunted hyperinflammatory responses. Induction of interferon responses within infected cells was rare and there was significant heterogeneity in the antiviral gene signatures, varying with the burden of infection in each cell. We also found that heavily infected secretory cells expressed abundant IL-6, a potential mediator of COVID-19 pathogenesis.

Introduction SARS-CoV-2, the virus responsible for COVID-19, primarily infects cells of the respiratory tract. The cellular tropism of SARS-CoV-2 may impact several aspects of the disease, including viral spread within and between hosts, mechanisms of immune control of infection or tissue pathology, and the therapeutic response to promising antivirals. Normal human tracheal bronchial epithelial (nHTBE) cells represent a diverse mix of ciliated epithelial cells, secretory cells, and basal cells that form a pseudostratified epithelium when cultured at the air-liquid interface, phenocopying the upper airway in humans [1,2]. Importantly, cells in this culture system also express endogenous levels of critical host factors including ACE2 and host proteases such as TMPRSS2 that are needed for SARS-CoV-2 viral entry [3][4][5][6][7]. This model also demonstrates key aspects of host antiviral epithelial immunity [8,9]. Recently, several studies using primary human lung cell cultures and respiratory cells isolated from SARS-CoV-2 infected patients have identified SARS-CoV-2 tropism for ciliated and secretory cells in the upper airway [10][11][12][13][14]. However, the heterogeneity of virus replication and induction of antiviral genes and proinflammatory cytokines within these cells is still unknown.
Remdesivir (GS-5734) has emerged as a promising direct antiviral therapy against SARS-CoV-2, with potent activity demonstrated against several coronaviruses [15,16]. A landmark clinical trial found that remdesivir treatment of hospitalized individuals with COVID-19 improved median recovery time [17], and this drug is now approved for COVID-19 under emergency use authorization by the U.S. Food and Drug Administration. Remdesivir is a prodrug that is metabolized in cells to the nucleotide analog remdesivir triphosphate, which interferes with coronavirus replication through delayed RNA chain termination [10,[18][19][20]. Recent studies have identified differential efficacy of remdesivir against SARS-CoV-2 in a range of cell culture systems linked to metabolism of the prodrug to the active form [10]. In addition to differential metabolism, other factors that may impact the variable efficacy of this drug in different cell types include differential drug uptake and heterogeneous permissibility of each cell type to viral entry and replication. While remdesivir clearly exhibits antiviral activity against SARS-CoV-2 in nHTBE cultures, it is not known if there are cell type-dependent differences in drug efficacy.
Following infection, coronaviruses are recognized by MDA5 and RIGI leading to the production of type I and III interferons (IFNs), which induce transcriptional programs that mobilize cellular antiviral defenses. Coronaviruses use several mechanisms to successfully evade detection resulting in rare and heterogeneous IFN production [21,22], similar to influenza virus infected cells [23][24][25][26]. During influenza A virus infection, we have previously identified interferon stimulated genes (ISGs) specifically induced in cells supporting high levels of virus replication and we have defined cell type-specific ISGs [27,28]. Additionally, we and others have found significant heterogeneity in antiviral responses across different cell types [28,29]. Cell type-specific responses and the degree of heterogeneity in antiviral responses can dictate the outcome of immune responses and infection.
Here, we use nHTBE cells infected with SARS-CoV-2 to demonstrate that remdesivir reduces viral replication uniformly in all susceptible cell types within the upper respiratory tract. Additionally, we demonstrate that TMPRSS2 is the primary host protease used for SARS-CoV-2 entry across cell types in the upper airway. Using single cell RNA sequencing, we further define SARS-CoV-2 tropism and the induction of antiviral and proinflammatory immune responses. We also uncover cellular heterogeneity in SARS-CoV-2 replication within the human airway epithelium and determine the relationship between productive viral infection and induction of IFN. Additionally, we define respiratory epithelial ISG expression in response to SARS-CoV-2 infection and uncover epithelial cell type-specific ISGs, and ISGs associated with either high levels of viral replication or protection from infection. Finally, we demonstrate that the proinflammatory cytokine IL-6 is primarily upregulated within heavily infected secretory cells. Together, these data help to define the early SARS-CoV-2 replication and antiviral landscape within the respiratory tract.

SARS-CoV-2 tropism and cell type-specific efficacy of remdesivir in differentiated tracheal bronchial epithelial cells
To determine how SARS-CoV-2 tropism varies within the upper airway we infected nHTBE cells differentiated at air-liquid interface with a SARS-CoV-2 infectious clone generated to express the fluorescent protein mNeonGreen (mNG) [30]. Using flow cytometry, we detected 12.5% of untreated nHTBE cells expressing mNG after icSARS-CoV-2-mNG infection (Fig 1A and 1B); treatment with remdesivir or IFNβ substantially reduced the percentage of mNG positive cells showing clear antiviral effect. We subsequently identified different cell types using the expression of several surface markers and SiR-tubulin ( Fig 1C) [31] and found that the ratio of different cell types in each culture remained similar during infection, regardless of anti-viral treatments ( Fig 1D). Among mNG + SARS-CoV-2 infected cells, we discovered a preponderance of ciliated (SiR-tubulin + CD271 -) and few infected secretory cells (CD66c + SiR-tubulin -CD271 -) (Fig 1E and 1F). We also identified a subset of infected cells lacking cell type-specific markers (termed "triple negative") and a population of cells expressing intermediate levels of CD271 and CD66c, but negative for SiRtubulin, termed "other". We did not identify mNG + basal cells (CD271 + SiR-tubulin -CD66c -). We also observed that ciliated cells expressed the highest mNG MFI, suggesting that this cell type is the most permissive for SARS-CoV-2 replication ( Fig 1G). These data are consistent with recent reports demonstrating ciliated cells as a dominant target cell of SARS-CoV-2 [10][11][12][13][14]32], and they suggest that basal cells are relatively resistant to infection. We next evaluated if the antiviral activity of remdesivir varies in different cell types depending on cellular uptake and metabolism [10], which could differ dramatically within the diverse population of cells in the upper respiratory tract. To determine how remdesivir impacts viral replication across all susceptible cells in the upper airway, we pretreated nHTBE cells with remdesivir at 0.1μM, a dose expected to inhibit 80% of SARS-CoV-2 replication [10] and infected with icSARS-CoV-2-mNG. As a positive control for antiviral activity, we also pretreated cells with IFNβ, which completely abrogated infection in nHTBE cells (Fig 1A and 1B), consistent with recent results demonstrating sensitivity of SARS-CoV-2 to interferon [33][34][35]. Remdesivir reduced both the number of SARS-CoV-2 infected cells by~5-fold and the expression level of mNG within infected cells (Fig 1B, 1F and 1G). Importantly, among infected cells, we found that remdesivir similarly inhibited replication regardless of cell type, reducing both the frequency of infection and the total mNG MFI (Fig 1F and 1G). These data suggest that cell type heterogeneity within the respiratory tract does not impact remdesivir efficacy.

Single cell RNA sequencing reveals heterogeneity in SARS-CoV-2 infection across differentiated tracheal bronchial epithelial cells
Tracheal-bronchial epithelial cells represent the first targets of SARS-CoV-2 and are key sentinels in the initiation of early antiviral responses. To further determine viral tropism and heterogeneity of replication and antiviral responses in the presence and absence of remdesivir, we infected nHTBE cells with SARS-CoV-2 and analyzed single cell transcriptomes via the 10X Genomics platform. We captured an average of 1729 cells per experimental condition and 1836 variable genes per cell. We observed significant heterogeneity in the level of virus replication supported per cell (Fig 2A, left). The frequency of cells with abundant viral RNA (defined as >0.1% of total reads per cell (S1A Fig) closely matched what we observed using flow cytometry (Fig 2A, right). Using this threshold for infection enables us to discriminate between cells that are productively infected from those with low-level viral RNA. Additionally, we confirmed viral replication kinetics with and without remdesivir by qRT-PCR, using a SARS-CoV-2 standard to calculate total viral genome copies (S1B Fig). Together these results show a distinct impact of remdesivir treatment on SARS-CoV-2 viral reads per cell, consistent with our flow cytometry results.

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia To determine SARS-CoV-2 tropism and the impact of remdesivir, we analyzed cells combined from all experimental conditions, including mock-infected and those at 24 or 48 hours post-infection (hpi), with or without remdesivir treatment, by uniform manifold approximation and projection (UMAP) [36,37]. While some clusters of cells in this analysis were defined by cell type (basal, ciliated, and ionocytes), the majority were undefined (Fig 2B and S1 Table). We further interrogated the transcriptional changes in these clusters by examining the gene ontology of the most significantly upregulated genes. This analysis confirmed the identify of ciliated and ionocytes cells and identified secretory cells as having a strong antiviral signaling response (S1C Fig). However, this analysis did not further illuminate the identity of the undefined cells. We observed several clusters that were characterized by a preponderance of infected cells ( Fig 2C); each of these clusters contained cells from different experimental conditions ( Fig 2D). For example heavily infected cells in clusters 2 and 7 are derived from cells infected at 24hpi, 24hpi with remdesivir and 48hpi. Consistent with our flow cytometry results, cells of a variety of phenotypes localized to each of these dominant infected cell clusters. Interestingly, the majority of cells treated with remdesivir (Fig 2A) clustered most closely with mock-infected cells. This clearly demonstrates the antiviral efficacy of remdesivir and suggests

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia that the presence of active SARS-CoV-2 replication is the strongest driver of transcriptomic phenotype in this cell population. These data suggest that SARS-CoV-2 infected airway epithelial cells adopt diverse transcriptional profiles dominated by both viral RNA level and cell type.

Single cell RNA sequencing identifies SARS-CoV-2 tropism in differentiated tracheal bronchial epithelial cells
Because UMAP clustering of combined samples was dominated by experimental condition, we next used t-squared neighbor embedded (t-SNE) clustering on individual samples to identify infected cell types and uninfected bystander cells (Fig 3A-3E and S2 Table). Ciliated cells were a major infected group at 24 and 48 hpi (Fig 3B and 3D), consistent with our flow cytometry data. We also observed heavily infected cells that did not cluster by cell type at 48 hpi (Fig 3D), possibly due to loss of cell markers during SARS-CoV-2 infection, as has also been observed with influenza infection [38]. To further evaluate the transcriptional changes in these clusters by examining the gene ontology of the most significantly upregulated genes. This analysis confirmed cell lineage for some clusters but was unable to further resolve the undefined cell types (S2A Fig). Similar to what others have observed [11,12,39], we found low level expression of ACE2 in our primary human epithelial cells (Fig 3A-3E), likely reflecting the limited sensitivity of scRNA sequencing. We therefore validated these findings using qRT-PCR and detected low ACE2 (30.4 ct) which increased after IFNβ treatment (26.07 ct,~25 fold) (S2B Fig), consistent with other studies [7,35]. Interestingly, we also observed augmented expression of the host protease TMPRSS2, implicated in SARS-CoV-2 entry [5], within ciliated cells, secretory cells, and ionocytes, but reduced or no expression in basal cells (Fig 3A-3E). To determine if TMPRSS2 is the primary protease utilized for SARS-CoV-2 entry in differentiated nHTBE cells we treated cells with the TMPRSS2 inhibitor camostat. Camostat has been shown to block pseudovirus infection in undifferentiated primary human lung [5] and SARS-CoV-2 infection in immortalized Calu3 cells [5,40]. Consistent with these results, blocking TMPRSS2 resulted in a near complete abrogation of infection in both ciliated and secretory cells (Fig 3F). These data demonstrate that TMPRSS2 is the dominant protease entry factor for cells of the upper airway and may help explain the relative resistance of basal cells to infection with SARS-CoV-2. nHTBE cells pre-treated with remdesivir had reduced numbers of infected cells (Fig 3C and  3E), similar to our flow cytometry data. Cell tropism of SARS-CoV-2 was not altered with remdesivir treatment, as infected cells were predominately ciliated cells and a few secretory cells at 24 hpi. Remdesivir pre-treatment had a stronger effect at 48 hpi as there were almost no infected cells detected.

Early antiviral immune response to SARS-CoV-2 in human airway epithelium
To evaluate the relationship between SARS-CoV-2 infection and antiviral immune responses we measured the induction of type I and III IFNs among individual cells in the human airway epithelium. To determine if there is a relationship between the levels of virus and the induction of IFN, we plotted any cell expressing IFNB, IFNL1, IFNL2, IFNL3 or IFNL4 by the level of SARS-CoV-2 infection. There was a strong positive correlation between levels of SARS-CoV-2 replication and induction of IFN (Fig 4A and 4B). Interestingly, we found cells with high levels of virus and low levels of IFN expression, potentially representing the successful evasion by SARS-CoV-2 of virus sensing pathways triggered by MDA5 or RIGI. At 48 hpi we observed two clusters of undefined cells, both with high SARS-CoV2 replication but divergent IFN responses: cluster 3 with high IFN induction and cluster 5 with reduced IFN. This disparate induction of IFN could indicate opposite outcomes of SARS-CoV-2 infection that drive

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia differential transcriptomic clustering. We identified very few (~0.5% at 24hpi) cells actively transcribing IFN, consistent with observations in other viral infections where a small number of IFN-producing cells is sufficient to establish the antiviral state [25,26]. Despite this, we detected robust induction of ISGs (Fig 4C). We evaluated ISGs with the greatest differences across sample condition and found several discrete patterns of expression. This analysis revealed ISGs that were only present in clusters of cells characterized with high levels of SARS-CoV-2 replication (Fig 4C, MICB and GEM). We also observed an ISG (MT1F) that was preferentially induced within ciliated cells (Fig 4C). At 48 hpi cells with the highest levels of viral RNA (clusters 3 and 5) had significantly lower levels of expression of several ISGs (DDIT, IFITM3, LY6E, TNFSF10), suggesting an impaired cellular antiviral immune response that enables unchecked SARS-CoV-2 replication (Fig 4C).
Notably with SARS-CoV-2 infection we observed induction of IL-6, a pro-inflammatory cytokine implicated in the pathogenesis of COVID-19 [41], only within infected cells at 24 and 48 hours post-infection. Heavily infected secretory cells expressed the highest levels of IL-6 at 24 hpi (Fig 4D). To determine the impact of early antiviral therapy on the hyperinflammatory response that drives COVID-19 pathology, we also analyzed infected cells pretreated with remdesivir. Remdesivir blocked the expression of IL-6 but not the induction of ISGs, suggesting that direct antiviral therapy can prevent induction of potentially pathogenic inflammation while preserving protective IFN responses.

Discussion
The SARS-CoV-2 pandemic continues to spread, causing significant worldwide morbidity and mortality. Understanding early events in virus-host interactions in relevant model systems will be critical for understanding the pathophysiology of the disease as well as for evaluating antiviral and immunomodulatory drugs. Here we demonstrate that SARS-CoV-2 has preferential tropism for ciliated cells and secretory cells in the upper airway. Within these primary cells SARS-CoV-2 is only sensed by a small number of infected cells leading to IFN activation. Despite the low number of cells expressing IFN, there is a robust antiviral ISG response established. Finally, we demonstrate that remdesivir is capable of restricting virus replication uniformly across all susceptible cells within the upper airway.
A recent study has demonstrated increased accumulation of remdesivir triphosphate in human airway epithelial cultures compared to Calu-3 and Vero cell cultures [18]. This study only addressed the accumulation and efficacy in readouts from total airway epithelial cell cultures. Here we are able to observe remdesivir efficacy in all of our individually defined cell populations using actual viral replication as the readout. In addition to blocking virus replication, remdesivir may also provide benefits from blocking induction of proinflammatory cytokines while maintaining antiviral signaling. We uncovered differential expression of the entry factor TMPRSS2 across multiple cell types with higher expression in ciliated cells and near absent expression in basal cells. The protease inhibitor camostat blocked SARS-CoV-2 infection across all susceptible cells, suggesting TMPRSS2 is the dominant protease that enables infection in the upper airway. Differential expression of this and other host entry factors may explain the resistance of basal cells to infection.

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia Coronaviruses employ diverse mechanisms to evade detection and induction of IFN. The exact mechanisms SARS-CoV-2 uses are unclear. However, SARS-CoV-2 clearly deploys successful strategies as induction of IFN within infected cells is rare (Fig 4A and 4B). While SARS-CoV-2 is proficient at evading IFN induction in most cells, it seems to be poor at countering effects of ISGs. Consistent with several recent reports, we demonstrate that SARS-CoV-2 is highly sensitive to IFN treatment (Fig 1) [33][34][35]. We have previously identified ISGs in epithelial cells associated with high levels of virus replication [28]. Congruently, cells with high levels of SARS-CoV-2 reads had significantly increased expression of some of these same ISGs including GEM (Fig 4D), RIPK2, PIM3, and ODC1 (S3 Table). We have previously demonstrated a number of cell type specific ISGs during influenza A virus infection in mice [28]. Interestingly, we also identified a cell type-specific ISG (MT1F) induced during SARS-CoV-2 infection, which was primarily expressed in ciliated epithelial cells. A better understanding of both the heterogeneity of ISG responses and cell-type specific induction will help to determine transcriptional combinations that are successful in blunting infection versus those that lead to establishment of infection, robust replication, and spread to new cells.
Together the results presented here demonstrate that ciliated airway epithelial cells are the predominant cell type initially infected by SARS-CoV-2 and reveal heterogeneity within infected cells and in antiviral responses. Importantly, remdesivir is capable of reducing viral replication in all infected cell types within this culture system. We also discovered that after viral spread at 48 hpi, IFN responses and TMPRSS2 expression are the major determinants of virus tropism in susceptible cells. Our results demonstrate that infected epithelial cells at the earliest stages can produce IL-6, which potentially contributes to the cascade of inflammatory disease that occurs later during infection. Defining virus-host interactions at the single cell level reveals broad tropism and successful evasion of virus sensing, key attributes that likely allowed the zoonotic emergence of SARS-CoV-2.

Materials and methods
Primary nHTBE cell culture nHTBE cells (Lonza Bioscience, CC-2540S) harvested from healthy individuals were expanded using Pneumacult-Ex Plus (Stem Cell Technologies) and seeded on 24 well hanging inserts, 0.4 μm pore (Corning, USA) and were maintained at ALI for 21 + days in a Pneumacult ALI growth medium (Stem Cell Technologies) at 37˚C/5% CO 2 in a humidified incubator. Cells were used 21 + days post ALI for experiments.

Drug treatments and infection with SARS-CoV-2
Basal layer of medium was replaced with infection media (DMEM with 2% FBS, 100 U/mL penicillin and streptomycin, 10 mM HEPES, 0.1 mM nonessential amino acids, 4 mM L-glutamate, and 1mM sodium pyruvate). Wells were left untreated or treated in the basal compartment with 100 ng of IFNβ or 0.1 μM remdesivir for 1 hr prior to infection, or with 50 μM of camostat mesylate for 2 hrs prior to infection. For infection, SARS-CoV-2 strain 2019-nCoV/ USA_WA1/2020 (provided by World Reference Center for Emerging Viruses and Arboviruses at the University of Texas Medical Branch) or icSARS-CoV-2-mNG, were added directly to the apical layer of cells for 1 hr at 37˚C with rocking every 10 min. The virus was removed from the apical compartment and cells were incubated at 37˚C until harvest.

Sample preparation for single-cell RNA sequencing
Cells were washed and trypsinized using ReagentPack Subculture Reagents (Lonza Bioscience). Single cell suspensions were washed with PBS + 5% FBS, resuspended in 200 μL PBS + 5% FBS, passed through a 70 μm filter, and placed on ice. Cells were counted using trypan blue exclusion on a Luna II (Logos Biosystems). The appropriate numbers of cells to achieve a targeted cell input of 5,000 cells per condition were used to generate the GEMs (Gel Bead-In Emulsion). The Chromium Next GEM Single Cell 3' Gel beads v3.1 kit (10X Genomics, Pleasanton, CA, USA) was used to create GEMs following manufacturer's instruction. All samples and reagents were prepared and loaded into the chip and ran in the Chromium Controller for GEM generation and barcoding. GEMs generated were used for cDNA synthesis and library preparation using the Chromium Single Cell 3' Library Kit v3.1 (10X Genomics) following the manufacturer's instruction.

Single-cell RNA sequencing analysis
The 10x Chromium Single Cell sequencing raw data of experimental samples were reference mapped using Cell Ranger Count (version 3.0.1) for alignment to a hybrid human+-SARS-CoV-2 viral genome reference (Homo_sapiens.GRCh38; MN985325). The hybrid genome reference, annotation file analysis scripts are available from our github site. The resulting raw count matrix for each experimental data set was imported into an R pipeline using Seurat 2.3.4 [42], where the filter criteria for empty droplets are minimum 1000 genes per cell, for genes that are presented in a minimum of three cells and for mitochondria genes percentage is no more than 30%. The default parameters were used for subsequent data normalization and scaling. The Seurat default method was used to select variable genes for PCA and clustering analysis, which the possible positive viral and mitochondria variable genes were excluded from final variable gene lists for all samples. The Seurat functions and utilities were used for subsequent cluster marker gene selection and results visualization. The ISG genes in the all marker list (FindAllMarkers function) of integrated analysis was ranked by their adjusted p value and ISGs within the top 30 that demonstrated disparate expression patterns are plotted with Seurat ViolinPlot function. To assess the relationship between all cells across all experimental conditions we used the FindClusters function with default parameters (except dim = 1:16 and resolution = 0.15) and the RunUMAP function with default parameters (dims = 1:16) of Seurat (version 3.2.2) was used for the integrated analysis for all five samples using only host genes. Gene ontology analysis was performed using Panther on genes in each cluster with an adjusted p value <0.01 and LogFC >0.5.

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia Overall cellular infection status (infected/uninfected) was determined by examining the distributions of percentage of viral counts within each cell on the log2 scale to magnify the differences at the low end. Thresholds for calling cells "infected" were set by calculating the percentage of summed viral counts within each cell over total counts, which is 0.1%. Cell types were defined using genes identified by Plasschaert et al [43]. A cutoff of at least 15 cell-typespecific genes was used to identify cell types for individual clusters.

PLOS PATHOGENS
Single cell analysis of SARS-CoV-2 tropism and antiviral immunity in primary airway epithelia