SARS-CoV-2 variants reveal features critical for replication in primary human cells

Since entering the human population, Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2; the causative agent of Coronavirus Disease 2019 [COVID-19]) has spread worldwide, causing >100 million infections and >2 million deaths. While large-scale sequencing efforts have identified numerous genetic variants in SARS-CoV-2 during its circulation, it remains largely unclear whether many of these changes impact adaptation, replication, or transmission of the virus. Here, we characterized 14 different low-passage replication-competent human SARS-CoV-2 isolates representing all major European clades observed during the first pandemic wave in early 2020. By integrating viral sequencing data from patient material, virus stocks, and passaging experiments, together with kinetic virus replication data from nonhuman Vero-CCL81 cells and primary differentiated human bronchial epithelial cells (BEpCs), we observed several SARS-CoV-2 features that associate with distinct phenotypes. Notably, naturally occurring variants in Orf3a (Q57H) and nsp2 (T85I) were associated with poor replication in Vero-CCL81 cells but not in BEpCs, while SARS-CoV-2 isolates expressing the Spike D614G variant generally exhibited enhanced replication abilities in BEpCs. Strikingly, low-passage Vero-derived stock preparation of 3 SARS-CoV-2 isolates selected for substitutions at positions 5/6 of E and were highly attenuated in BEpCs, revealing a key cell-specific function to this region. Rare isolate-specific deletions were also observed in the Spike furin cleavage site during Vero-CCL81 passage, but these were rapidly selected against in BEpCs, underscoring the importance of this site for SARS-CoV-2 replication in primary human cells. Overall, our study uncovers sequence features in SARS-CoV-2 variants that determine cell-specific replication and highlights the need to monitor SARS-CoV-2 stocks carefully when phenotyping newly emerging variants or potential variants of concern.


Introduction
Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), a novel betacoronavirus belonging to the Coronaviridae family, appears to have first entered humans in late 2019 in the Hubei Province of China [1]. Since then, it has spread worldwide in the human population, predominantly causing a mild to severe respiratory disease, termed "Coronavirus Disease 2019 ." Currently, pandemic SARS-CoV-2 has led to more than 100 million laboratory-confirmed COVID-19 cases globally, resulting in more than 2 million deaths [2]. Coronaviruses are enveloped viruses with large positive-sense RNA genomes of approximately 30 kb [3]. The SARS-CoV-2 Spike (S) glycoprotein resides on the surface of virions and mediates viral entry into the host cell by binding to cellular ACE2 [4,5] and triggering viralhost membrane fusion [3,6]. This fusion function of S is dependent on its cleavage by host cell proteases, which occurs either following attachment of virions to the host cell membrane or during virion maturation and egress [3,7,8]. Strikingly, SARS-CoV-2 S harbors a furin cleavage site [9], and pre-activation "priming" of S by furin-like proteases increases fusion-dependent entry efficiency of the virus at the plasma membrane following further cleavage by host TMPRSS2 [10][11][12][13]. Other SARS-CoV-2 virion components include the Membrane (M) protein that confers shape and support to the virus particle and interacts with both S and the Nucleoprotein (N; which coats the viral RNA genome [14,15]) and the Envelope (E) protein, a small membrane protein that both drives virion assembly and budding [16] and has ion channel activity linked to viral pathogenesis [17].
Following S-mediated entry, translation of the viral genomic RNA results in expression of several nonstructural proteins (nsp1 to 16), many of which form the essential replicase-transcriptase complex or fulfill additional important functions in the virus life cycle, such as evasion of host innate immunity [3]. Notably, the replication complex of coronaviruses includes proofreading activity, leading to greater genome stability as compared to other RNA viruses that typically lack this function [18]. Nevertheless, several mutations in the SARS-CoV-2 genome have already been reported in human-circulating viruses worldwide, and some exhibit increasing prevalence suggestive of positive selection or fitness advantage in the new human host [19][20][21]. An important example of this is the D614G substitution in S, which rapidly dominated SARS-CoV-2 sequences following its emergence, and is functionally linked to increased virus infectivity and replication capacity [22][23][24][25], potentially aiding virus spread throughout the population. Other recently acquired mutations in the SARS-CoV-2 genome also appear to confer an advantageous replication and/or transmission phenotype in humans and may impact virus antigenicity and the effectiveness of deployed vaccines [21].
In this study, we characterized 14 SARS-CoV-2 isolates from patient material collected during the first pandemic wave in Switzerland between March and May 2020. Phylogenetically, the isolates were representative of the different virus clades circulating in Europe in early 2020 [26]. While the individual isolates each carried distinct genetic variants that were already detected in the original patient material, rare additional mutations (notably in S and E) were observed upon low-passage virus stock preparation in Vero-CCL81 cells. Comparative replication and passage analysis of each SARS-CoV-2 isolate in Vero-CCL81 cells and differentiated primary human bronchial epithelial cells (BEpCs) revealed the critical importance of the S furin cleavage site and E positions 5/6 as specific determinants of SARS-CoV-2 infectivity in human respiratory cells. Furthermore, we observed that SARS-CoV-2 isolates expressing S G614 generally appeared to replicate more efficiently than isolates harboring S D614. Notably, SARS-CoV-2 isolates with naturally occurring Orf3a (Q57H) and nsp2 (T85I) substitutions exhibited a cell-specific replication phenotype, with efficient propagation restricted to primary human BEpCs. Overall, our comparative functional and sequence characterization of a broad range of virus isolates describes both important conserved residues, as well as naturally occurring substitutions, which contribute to efficient SARS-CoV-2 replication in human respiratory cells.

Isolation and sequence analysis of SARS-CoV-2 from patient diagnostic samples in Switzerland, March to May 2020
Following identification of the first SARS-CoV-2-positive patient in Switzerland at the end of February 2020, the initial pandemic wave was characterized by a peak of >1,000 new confirmed cases per day by mid-March and a slow decline to a relatively stable approximately 50 new cases per day by mid to late May (Fig 1A). From a biobank of over 1,650 SARS-CoV-2 reverse transcription PCR (RT-PCR)-positive patient samples collected for diagnostic purposes throughout this period, we attempted virus isolations from 67 independent samples. Isolation success could readily be stratified by RT-PCR cycle threshold (Ct) value: We were able to isolate 62% of SARS-CoV-2 isolates with initial diagnostic Ct values <25 (n = 21), while SARS-CoV-2 was only isolated from a single sample with diagnostic Ct values >25 (n = 46) (Fig 1B, S1 Table). These isolation data are broadly in line with SARS-CoV-2 culture success rates and correlations with Ct values described by others [24,[27][28][29] and indicate that presence of cell culture-infectious virus can be predicted from initial diagnostic RT-PCR results.
Fourteen SARS-CoV-2 primary isolates (IMV1-14) were cultured from samples obtained between March 10, 2020 and May 4, 2020 (highlighted on Fig 1A). Together with BavPat1, a SARS-CoV-2 isolate collected on January 29, 2020 in Munich, Germany [30], all viruses were minimally propagated on Vero-CCL81 cells to obtain working stocks (passage 3 (P3) for IMV1-14; passage 2 (P2) for BavPat1). Original Swiss patient material, as well as both P2 and P3 stocks from all isolates, was subjected to next-generation sequencing (NGS) to confirm the absence of additional pathogens and to provide full-length SARS-CoV-2 genomic sequences (S2 and S3 Tables). Phylogenetic analysis based on the P2 sequences revealed that the isolated viruses were representative of viruses found across Europe during the first half of 2020 [26] (S1 Fig). As expected, all 14 SARS-CoV-2 isolates harbored several amino acid substitutions when compared to the reference strain SARS-CoV-2/Wuhan-Hu-1, derived from a patient in Wuhan, China in December 2019 [31]. We noted that during working stock preparation (i.e., a total of 3 passages from original patient material in Vero-CCL81 cells), some of the virus isolates acquired additional amino acid substitutions or deletions that were not detected in the original patient material nor in early (P2) passages (S2 and S3 Tables). Most of these passageinduced amino acid substitutions were observed in S, E, or nsps. However, there was no clear consensus change across all isolates that were passaged similarly. The most notable changes were in E, where 6 out of 14 isolates exhibited amino acid substitutions following working stock growth in Vero-CCL81 cells (S3 Table). Furthermore, passage-induced deletion of the furin cleavage site in S was only observed in 2 out of the 14 isolates (IMV1: Δ677-688, 17%; IMV14: Δ679-685, 81%) and was fully retained in the other 12 isolates. It is possible that the low number of passages that we performed to generate working stocks limited the selection of S furin cleavage site deletions that have been otherwise readily detected by others [32][33][34]. We conclude that our set of 14 isolates (IMV1-14) is representative of viruses circulating in Europe in early 2020. Together with their lack of consensus-adaptation to Vero-CCL81 cells during stock preparation, these isolates should therefore be suitable to functionally characterize potential features of early SARS-CoV-2 human adaptation and replication.

Characterization of SARS-CoV-2 isolates in Vero-CCL81 cells and primary human bronchial epithelial cells
Using Vero-CCL81 cells, we assessed the plaque phenotypes and growth kinetics of all 14 SARS-CoV-2 isolates, as well as BavPat1. All viruses plaqued well on this substrate, although there was some heterogeneity in plaque size between isolates: For example, as compared to BavPat1, it was striking that IMV4, IMV8, and IMV11 produced smaller plaques (Fig 1C and  1D). Differences between isolates were also observed during multi-cycle growth analysis. While most viruses grew to titers around 10 6 or 10 7 plaque-forming unit (PFU)/mL over 72 h, some viruses were clearly attenuated: Final titers of IMV2, IMV8, IMV9, and IMV14 were 10to 100-fold lower than the other isolates, and kinetic analysis showed that IMV1, IMV3, and IMV13 grew slower than other isolates (Fig 1E and 1F). Such phenotypic differences were not specific to Vero-CCL81 cells, as we also assessed the replication of a subset of isolates in Vero-E6 cells and noted similar viral growth patterns (S2 Fig).
We next assessed SARS-CoV-2 isolate replication in a previously described primary human BEpC model (Fig 2A, [35]). BEpCs were grown at air-liquid interface (ALI) for a minimum of 4 weeks, and their differentiation into a pseudostratified respiratory epithelium was validated by measuring increased transepithelial electrical resistance (TEER) and epithelium-specific cell and tight junction markers, such as β-tubulin and zona occludens protein 1 (ZO-1) (S3A and S3B  Fig). BEpCs were infected from the apical side with each SARS-CoV-2 isolate, and viral titers in both the apical and basolateral compartments were monitored over 72 h. Strikingly, there were clear differences in viral replication kinetics between the individual isolates in the apical washes of the infected epithelium (Fig 2B and 2C). BavPat1, as well as 5 other isolates (IMV2, IMV7, IMV9, IMV10, and IMV12), grew to high titers (>10 7 PFU/mL) within 72 h. In stark contrast, 2 isolates (IMV6 and IMV4) were unable to replicate in BEpCs at all, while IMV14 was strongly attenuated in this tissue substrate, only reaching titers below 10 4 PFU/mL. IMV11 also exhibited a notably attenuated growth phenotype, yielding 100-fold lower titers compared to BavPat1, while IMV1, IMV3, IMV5, IMV8, and IMV13 exhibited intermediate growth phenotypes (approximately 10-fold lower than BavPat1) (Fig 2B and 2C). SARS-CoV-2 titers in the basolateral compartment were minimal as compared to titers in the apical compartment (S3C Fig), and virus was only ever detectable after 72 h of infection, suggesting damage of the bronchial epithelium in response to SARS-CoV-2. We used NGS to fully sequence all viral isolates from BEpC apical compartment samples taken following 72 h of replication, but did not observe any mutations leading to amino acid substitutions as compared to the input Vero-derived P3 stocks used for initial infection. Furthermore, to validate our observations, we generated differentiated BEpCs from 2 additional human donors and used these new primary cultures to assess the replication capabilities of a subset of virus isolates. While there was some donor-to-donor variation in peak viral titers, most isolates grew very well, but IMV6 and IMV14 were clearly attenuated in all donors, and IMV4 was attenuated in 2 out of 3 donors (S4 Fig). [www.bag.admin.ch]). Circles depict the day of patient sample acquisition of the 14 SARS-CoV-2 isolates generated in this study. (B) SARS-CoV-2 isolation success rate stratified by initial diagnostic Ct value of patient material. A total of 67 PCR-positive patient samples were subjected to virus isolation attempts on Vero-CCL81 cells. Positive virus cultures exhibited cytopathic effect and low SARS-CoV-2-specific Ct values (<22) after 4-6 days of culture. Numbers at the top of bars represent isolation success/attempts. Numbers on x-axis represent Ct values. (C) Plaque phenotypes of all SARS-CoV-2 isolate working stocks (BavPat1, IMV1-14) on Vero-CCL81 cells. Monolayers were infected with similar PFUs of each isolate and cultured with an agar overlay for 3 days prior to fixing and staining with crystal violet. Shown is a representative image from each isolate used. (D) Quantification of plaque sizes from (C). Individual data points are shown, with means and standard deviations represented. (E) Vero-CCL81 cells were infected with the indicated SARS-CoV-2 isolates at an MOI of 0.01 PFU/cell. Supernatants were harvested at several times postinfection, and SARS-CoV-2 titers were determined by plaque assay. Data represent means and standard deviations from 3 independent experiments, with each experiment performed using technical duplicates. The dotted line crossing the y-axis at 10 2 PFU/mL indicates the assay LoD. (F) AUC values for the virus replication data shown in (E), normalized to inoculum titer of the respective isolate. Data represent means and standard deviations from the 3 independent experiments. A 1-way ANOVA was performed on log2-transformed data to test for significant differences between the individual isolates (p < 0.0001). To test for statistical significance against BavPat1, an unpaired t test was used on log2-transformed data ( � p < 0.05; �� p < 0.01; ��� p < 0.001). For underlying data, see S1 Data. AUC, area under the curve; Ct, cycle threshold; LoD, limit of detection; MOI, multiplicity of infection; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2. https://doi.org/10.1371/journal.pbio.3001006.g001

PLOS BIOLOGY
SARS-CoV-2 variants reveal features critical for replication in primary human cells were infected with 6,000 PFU of each SARS-CoV-2 isolate from the apical side. At several times postinfection, apical washes were harvested, and virus titers were determined by plaque assay. Data represent means and standard deviations from 3 in dependent experiments, with each experiment performed using 1 individual transwell. The dotted line crossing the y-axis at 10 2 PFU/mL indicates the assay LoD. (C) AUC values for the virus replication data shown in (B), normalized to inoculum titer of the respective isolate. Data represent means and standard deviations from the 3 independent experiments. A 1-way ANOVA was performed on log2-transformed data to test for significant differences between the individual isolates (p < 0.0001). To test for statistical significance against BavPat1, an unpaired t test was used on log2-transformed data ( � p < 0.05; ��� p < 0.001). (D) Indirect immunofluorescent staining of SARS-CoV-2-infected BEpCs. Differentiated BEpCs were infected from the apical side with 6,000 PFU of the indicated SARS-CoV-2 isolate. At 72 We next performed indirect immunofluorescent staining in the differentiated human BEpCs to determine the cell tropism of a subset of viral isolates that represented the various growth phenotypes we had previously observed. Using antibodies against β-tubulin, mucin 5AC, uteroglobin, and P63, we identified ciliated, goblet, club, and basal cells, respectively. Viral infection was monitored using an anti-double-stranded (ds) RNA antibody. Despite their reproducible differences in replication capabilities, we found that all isolates displayed a similar cell tropism in this BEpC model of the human respiratory epithelium that was independent of donor source: Similar to previous reports [36,37], we observed SARS-CoV-2 infection of secretory (goblet and club) cells, as well as ciliated cells, although we note that SARS-CoV-2-infected β-tubulin-positive cells often displayed reduced cilia integrity, which has recently been reported by others [38,39] (S5-S10 Figs). Overall, these SARS-CoV-2 functional data in Vero-CCL81 cells, Vero-E6 cells, and 3 independent BEpC donors reveal striking differences in replication properties between viral isolates despite a similar cell tropism in primary differentiated human respiratory cells.

Host cell-specific replication efficiency of SARS-CoV-2 isolates
To understand if the spectrum of SARS-CoV-2 isolate replication kinetics was host cell specific, we compared the relative replication of individual isolates in Vero-CCL81 cells and BEpCs ( Fig  3A). While some virus isolates replicated very well in both cell types (BavPat1, IMV10, and IMV12), other isolates displayed striking host-specific phenotypes. Firstly, IMV2, IMV8, and IMV9 replicated relatively well in BEpCs, but were highly attenuated in Vero-CCL81 cells ( Fig  3A). These 3 isolates were unique in this phenotype, and each harbored combined amino acid substitutions in Orf3a (Q57H) and nsp2 (T85I) that were not identified in any other isolate (S3 Table), suggestive of an association between these specific changes and the observed cell-specific replication. In contrast, IMV4, IMV6, and IMV11 exhibited a high replication capacity in Vero-CCL81 cells, but were strongly attenuated in BEpCs, with IMV4 and IMV6 being particularly attenuated in most donors (Fig 3A, S4E Fig). Notably, these 3 isolates each harbored amino acid substitutions at positions 5 (V5G/A) or 6 (S6W) of E (S3 Table). We also noted that IMV14 was highly attenuated in BEpCs and moderately attenuated in Vero-CCL81 cells, which could be due to >80% of the virus population containing a deletion of the furin cleavage site in S (Δ679 to 685). As observed by others [24], SARS-CoV-2 isolates such as BavPat1, IMV5, IMV7, IMV10, IMV12, and IMV13 that express S G614 replicated more efficiently in BEpCs than IMV1 and IMV3 that harbor S D614 (Fig 3A, S3 Table). These comparative replication data across 14 different SARS-CoV-2 isolates suggest associations between specific amino acid residues and distinct substrate phenotypes and imply that the functional viral traits underlying these residues are likely to have important roles for virus replication in human respiratory cells.

Prevalence of functional SARS-CoV-2 amino acid substitutions worldwide and in laboratory stocks
We investigated the prevalence of the phenotype-associated SARS-CoV-2 isolate variants in Orf3a, nsp2, E, and S in our different Vero-CCL81 passages, original patient material, and h postinfection, cells were fixed, permeabilized, and stained for the presence of infected cells (dsRNA; magenta) and ciliated cells (βtubulin; cyan), goblet cells (Muc5AC; cyan), club cells (uteroglobin; cyan), or basal cells (P63; yellow). Nuclei were stained with DAPI (blue). Z-stacks were transformed into maximum projection images. Scale bar represents 9.6 μm. Arrows indicate co-localization. Data show a subset of representative images from 2 independent experiments performed in different BEpC donors (see also S5-S10 Figs). For underlying data, see S1 Data. ALI, air-liquid interface; AUC, area under the curve; BEpC, bronchial epithelial cell; LoD, limit of detection; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2.
https://doi.org/10.1371/journal.pbio.3001006.g002 For WW prevalence, published consensus sequences were used, while for patient samples and passages, the available NGS data were analyzed, and prevalence was determined based on a frequency of >15%. (E and F) Variant frequency of selected SARS-CoV-2 E and S variants within individual patient samples and Vero-CCL81 passages for the indicated isolate was determined by NGS. For worldwide. Orf3a (Q57H) and nsp2 (T85I) variants were found in 20% of our patient samples, a similar prevalence to that found worldwide for the individual variants (Fig 3B). Notably, this proportion did not change upon passaging of isolates in Vero-CCL81 cells during stock preparation. In contrast, E protein variants at positions 5/6 and furin cleavage site deletions in S were all undetectable in worldwide and patient samples, but clearly became more prevalent among our SARS-CoV-2 isolates upon passaging in Vero-CCL81 cells (Fig 3C and 3D). This passaging effect was further exemplified when reanalyzing NGS data from each individual isolate and its parental patient material: The frequency of E S6W, V5G, or V5A increased substantially during stock preparation of IMV4, IMV6, and IMV11, respectively (Fig 3E), and similar increases in furin cleavage site deletion variants were observed during stock preparation of IMV1 and IMV14 (Fig 3F). These data indicate that, unlike the Orf3a (Q57H) and nsp2 (T85I) variants that exhibit biased replication to BEpCs, the E (V5G/A, S6W) and S furin cleavage site deletion variants that exhibit biased replication to Vero-CCL81 cells are positively selected for during passaging in Vero-CCL81 cells.

Comparative passaging reveals importance of the spike furin cleavage site for SARS-CoV-2 replication in primary human bronchial epithelial cells
The furin cleavage site in SARS-CoV-2 S is highly conserved in patient samples worldwide (Fig 3D) and is clearly important for virus pathogenesis [40]. However, a number of studies, including our own here, have reported that deletions in the furin cleavage site can be acquired during passage in Vero cell substrates [32,33]. To experimentally assess the importance of the SARS-CoV-2 S furin cleavage site in primary human respiratory cells, we took advantage of 2 isogenic IMV14 isolate stocks that exhibited differing frequencies of the S furin cleavage site deletion (P2, 23.5%; and P3, 81.1%, Fig 4A). Comparative passaging for 4 independent replicates in Vero-CCL81 cells and BEpCs revealed that IMV14 P2 (23.5% S furin cleavage site deletion) replicated in both cell systems (Fig 4B and 4D). Furthermore, sequencing of each passage supernatant revealed that, in Vero-CCL81 cells, IMV14 P2 generally increased its frequency of S furin cleavage site deletion from approximately 20% to 50% to 80% (3 out of 4 replicates within 2 additional passages; Fig 4C). In contrast, there was a rapid loss of this deletion variant during passage in BEpCs, with 100% of recovered sequences for all 4 replicates harboring an intact S furin cleavage site after 2 additional passages (Fig 4E). Broadly, similar results were obtained when using IMV14 P3 (81.1% S furin cleavage site deletion): The virus replicated well in Vero-CCL81 cells and retained the S furin cleavage site deletion at approximately 80% frequency (Fig 4F and 4G). However, IMV14 P3 was highly attenuated in BEpCs, and virus was only recovered in sufficient quantities from 1 of 4 replicates (Fig 4H). Strikingly, sequencing of this BEpC-recovered virus revealed a rapid loss of the S furin cleavage site deletion (Fig 4I). Together with work from others, these data provide strong experimental evidence that the S furin cleavage site is an essential viral trait required, and selected, for efficient SARS-CoV-2 replication in primary human respiratory cells [40,41].

PLOS BIOLOGY
SARS-CoV-2 variants reveal features critical for replication in primary human cells Firstly, using fully replication-competent SARS-CoV-2 isolates in primary human BEpCs, we observe that viruses expressing the naturally occurring G614 substitution in S generally replicate more efficiently than otherwise similar viruses expressing the parental D614 version of S. These observations are in line with previous work using virus-like particle or vesicular stomatitis virus (VSV) pseudotype systems, where S G614 pseudotypes were found to be more infectious than S D614 pseudotypes [22,23]. Furthermore, recent SARS-CoV-2 reverse genetics studies have also linked the single D614G substitution to increased infectivity and virus replication in primary human cells as well as in vivo models [24,25]. Mechanistically, the D614G substitution appears to permit a more open ACE2-binding conformation of the S trimer, which may lead to increased virus fusion with host cell membranes [22]. Together with our results from naturally occurring SARS-CoV-2 isolates, these common findings using different experimental systems provide a clear basis to understand the rapid emergence and populationwide spread of the S D614G substitution.
Secondly, we uncovered an association between naturally occurring Orf3a (Q57H) and nsp2 (T85I) variants, and a cell-specific replication phenotype: SARS-CoV-2 isolates unique in harboring these variants replicated poorly in Vero-CCL81 cells, while maintaining efficient replicative capacity in primary human BEpCs. These variants are found in approximately 20% of SARS-CoV-2 sequences worldwide, but their functional consequences are unclear. The functions of SARS-CoV-2 nsp2 have not been fully elucidated, although it has been described to interact with host proteins involved in endomembrane compartments, vesicle trafficking, and translation [42]. Orf3a is an ion channel with proapoptotic and pro-inflammatory functions [42][43][44]. Recent structural studies have revealed that Q57 forms the major hydrophilic constriction in the SARS-CoV-2 Orf3a pore; however, the Q57H substitution does not appear to influence Orf3a channel activities [44]. Nevertheless, the Orf3a Q57H variant was acquired shortly after introduction into humans and is not found in early SARS-CoV-2 sequences or in related bat coronaviruses [44], potentially suggesting positive selection of this residue similar to S D614G. Notably, the mutation leading to the Orf3a Q57H variant also changes the sequence of 2 other putative (and poorly characterized) overlapping reading frame products, Orf3c and Orf3d [45,46]. Further studies, including reverse genetics experiments, will be required to determine whether these variants are functionally linked to one another and how each variant may mechanistically determine this phenotype. While Vero-CCL81 cells are clearly not a physiological substrate for SARS-CoV-2, it may be that the cell-specific replication capacities that we observe with isolates expressing Orf3a (Q57H) and nsp2 (T85I) variants also occurs in other, non-respiratory, cell types of the human body, potentially impacting disease pathogenesis.
Thirdly, we observed that SARS-CoV-2 stock preparation in Vero-CCL81 cells led many of the isolates (6 out of 14) to acquire passage-derived amino acid substitutions in the E protein, which were not observed in patient material (either from our samples or worldwide). Strikingly, substitutions at E positions 5/6 (V5G/A or S6W) were highly attenuated in primary human BEpCs but replicated well, and were potentially selected for, in Vero-CCL81 cells. Previous work using SARS-CoV-1 has shown that E deletion viruses are only mildly attenuated in independent replicates. (F-I) IMV14 P3 stock (81.1% S furin cleavage site deletion frequency) was passaged twice on Vero-CCL81 cells or BEpCs. At each passage, virus titers were determined by plaque assay (F and H), and frequency of the furin cleavage site deletion was determined by NGS (G and I). Data represent the individual results from 4 independent replicates. For E and I, a sequence could only be determined from 1 replicate, where limited SARS-CoV-2 replication occurred. In panels B, D, F, and H, the dotted lines crossing the y-axes at 10 2 PFU/mL indicate the assay LoD. For underlying data, see S1 Data. BEpC, bronchial epithelial cell; LoD, limit of detection; NGS, next-generation sequencing; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2. https://doi.org/10.1371/journal.pbio.3001006.g004

SARS-CoV-2 variants reveal features critical for replication in primary human cells
Vero-E6 cells as compared to their attenuation in human cell lines or in in vivo models [16], and such viruses have been considered as the basis for live-attenuated vaccine concepts [47][48][49][50]. These data, together with the observations presented here, suggest that SARS-CoV-2 E has important species-or cell type-specific functions and further indicate that residues at positions 5 and 6 play key roles in some critical E activities. It will be interesting in future studies to determine whether such activities include slowing of the cell secretory pathway by E [51], and/ or the activation of host inflammasomes by E [52], as well as how the requirement for these functions differs between species or cell types.
Finally, as also observed by others [32][33][34]53], we found that in-frame deletions of the furin cleavage site in S could be selected for and enriched in SARS-CoV-2 isolates during passage in Vero-CCL81 cells. This was apparent even for very low-passage stocks as described here (passage 3 from original patient material), but importantly was not a universal feature of all SARS--CoV-2 isolates: Deletions of the S furin cleavage site were only detected in 2 out of 14 of our primary isolates. The furin cleavage site permits "priming" of S (cleavage at S1/S2 boundary) by furin-like proteases [11], thereby enabling the virus to fuse at the cellular plasma membrane if TMPRSS2 is present to further cleave S at the S2' site [12,13]. In the absence of furin cleavage, "priming" and subsequent cleavage of S at the S2' site occurs in endosomes by cathepsin proteases. Vero cells appear to lack expression of TMPRSS2 [13,33], and thus maintaining the S furin cleavage site is unlikely to be advantageous to SARS-CoV-2 in these cells. This provides a plausible mechanistic basis for SARS-CoV-2 adaptation to Vero cells through deletions in the S furin cleavage site. Such findings, together with our additional observations relating to functional E substitutions during SARS-CoV-2 growth in Vero-CCL81 cells, underscore the need for care during stock preparation, particularly where viruses are to be compared with one another for phenotypes. Thus, it may be that other cell lines, including human Calu-3 [53], are more suitable for SARS-CoV-2 propagation and should be thoroughly characterized in the future as a possible substrate to limit cell culture adaptation in SARS-CoV-2 genes such as S and E. Nevertheless, we found that in-frame deletions in the S furin cleavage site attenuated SARS-CoV-2, most notably in primary human BEpCs. However, passaging of SARS-CoV-2 preparations containing the S furin cleavage site deletions in primary human BEpCs provided a very high selection pressure to enrich only for viruses that retained the full S furin cleavage site. Thus, our study of SARS-CoV-2 cell adaptations revealed the critical nature of this S furin cleavage site for SARS-CoV-2 replication in primary human respiratory cells, a finding similar to that recently described by others using alternative methods [40,41].
In sum, our work identifies both key features in the SARS-CoV-2 genome that are essential for its growth in human respiratory cells, as well as variants derived during human circulation that determine efficient replication.

SARS-CoV-2 full genome sequencing
SARS-CoV-2 genome sequencing from P2 and P3 stocks was performed using an NGS approach as described previously [57]. Briefly, samples were diluted in NucliSENS EasyMAG lysis buffer (BioMérieux, Craponne, France; ratio 1:4) and filtered using a 0.45-μm PES filter (TPP, Trasadingen, Switzerland). Total nucleic acids were extracted using the NucliSENS EasyMAG system, followed by reverse transcription with random hexamers and second strand synthesis. Sequencing libraries were constructed using the NexteraXT protocol (Illumina, San Diego, California, USA) and sequenced on an Illumina MiSeq for 1 × 151 cycles using version 3 chemistry.
SARS-CoV-2 genome sequencing from original patient material was performed using a previously described targeted, tiled-amplicon, NGS approach [58,59]. Briefly, total nucleic acids were extracted followed by reverse transcription with random hexamers and oligo-dT priming (ratio 3:1) using SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific, USA). The generated cDNA was used as input for 14 overlapping PCR reactions (ca. 2.5 kb each) spanning the entire viral genome using Platinum SuperFi DNA Polymerase (Thermo Fisher Scientific). Amplicons were pooled per patient before NexteraXT library preparation and sequencing on an Illumina MiSeq for 1 × 151 cycles.

SARS-CoV-2 replication assays
For Vero-CCL81 and Vero-E6 experiments, 1.2 × 10 4 cells were seeded overnight in 96-well plates. Cells were then infected with SARS-CoV-2 at an MOI of 0.01 PFU/cell in PBS supplemented with 0.3% BSA, 1 mM Ca 2+ /Mg 2+ , 100 U/mL penicillin, and 100 μg/mL streptomycin. After 1-h inoculation, cells were washed once in PBS, and the medium was replaced with DMEM supplemented with 100 U/mL penicillin, 100 μg/mL streptomycin, 0.3% BSA, 20 mM HEPES, 0.1% FCS, and 0.5 μg/mL TPCK-treated trypsin. Samples were taken at the indicated time points and stored at −80˚C prior to titer determination by plaque assay [35]. Plaque size quantification was performed using ImageJ analysis software. For BEpC experiments, cells grown on 6.5-mm transwell filter inserts were infected from the apical side with 6,000 PFU of the respective SARS-CoV-2 isolate in PBS supplemented with 0.3% BSA, 1 mM Ca 2+ /Mg 2+ , 100 U/mL penicillin, and 100 μg/mL streptomycin. After 1 h, the inoculum was removed, and cells were incubated at 37˚C. At selected time points, virus was harvested from the apical side by applying 80 μL of PBS for 15 min at 37˚C before removal and storage at −80˚C prior to titer determination by plaque assay [35]. Basolateral samples were also collected and frozen at −80˚C.

SARS-CoV-2 passaging assays
For Vero-CCL81 experiments, approximately 2 × 10 5 cells were infected with the indicated SARS-CoV-2 isolate stock at an MOI of 0.01 PFU/cell as described above. In parallel, BEpCs grown on a 12-mm filter insert were infected from the apical side with 2.2 × 10 4 PFU/well. After 72 h at 37˚C, virus was harvested by either directly collecting the supernatant from infected Vero-CCL81 cells or by performing apical washes of BEpCs with PBS. One-third of the harvested supernatant (passage 1; P1) was used to similarly infect fresh Vero-CCL81 cells or BEpCs for a subsequent passage (P2), while the remaining supernatant was used for full genome sequencing (see above) and SARS-CoV-2 titration [35].

Statistical analyses
Statistical significance was determined using 1-way ANOVA and unpaired t test on log2-transformed data. , normalized to inoculum titer of the respective isolate. Data represent means and standard deviations from the 3 independent experiments. A 1-way ANOVA was performed on log2-transformed data to test for significant differences between the individual isolates (p < 0.0001). To test for statistical significance against BavPat1, an unpaired t test was used on log2-transformed data ( �� p < 0.005; ��� p < 0.0005). For underlying data, see S1 Data. AUC, area under the curve; LoD, limit of detection; MOI, multiplicity of infection; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2. nuclei (DAPI; blue). Scale bar represents 25 μm. (C) Differentiated BEpCs from donor 1 were infected from the apical side with 6,000 PFU of each SARS-CoV-2 isolate (see Fig 2B). At the indicated times postinfection, basolateral samples were harvested, and virus titers were determined by plaque assay. Data represent means and standard deviations from 2 independent replicates. The dotted line crossing the y-axis at 10 2 PFU/mL indicates the assay LoD. For underlying data, see S1 Data. ALI, air-liquid interface; BEpC, bronchial epithelial cell; LoD, limit of detection; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2; TEER, transepithelial electrical resistance; ZO-1, zona occludens protein 1. , normalized to inoculum titer of the respective isolate. Data represent means and standard deviations from the 3 independent experiments. A 1-way ANOVA was performed on log2-transformed data to test for significant differences between the individual isolates (n.s. for donor 2; p < 0.0001 for donor 3). To test for statistical significance against BavPat1 (D), an unpaired t test was used on log2-transformed data ( � p < 0.05; �� p < 0.01). Of note, as BavPat1 exhibited an unexplained attenuated phenotype in BEpCs from donor 2 (B), an unpaired t test was used on log2-transformed data to test for statistical significance against IMV3 ( � p < 0.05; �� p < 0.01). (E) Boxplot representation of the AUC values for the virus replication data from the 3 independent donors (data from B and D; and Fig 2C), normalized to inoculum titer of the respective isolate and donor. Shown is median and standard deviation. The red triangles represent data obtained from donor 1. For underlying data, see S1 Data. AUC, area under the curve; BEpC, bronchial epithelial cell; LoD, limit of detection; PFU, plaque-forming unit; SARS-CoV-2, Severe Acute Respiratory Syndrome Coronavirus 2.