COVID-19 global cases have climbed to more than 33 million, with over a million total deaths, as of September, 2020. Real-time massive SARS-CoV-2 whole genome sequencing is key to tracking chains of transmission and estimating the origin of disease outbreaks. Yet no methods have simultaneously achieved high precision, simple workflow, and low cost. We developed a high-precision, cost-efficient SARS-CoV-2 whole genome sequencing platform for COVID-19 genomic surveillance, CorvGenSurv (Coronavirus Genomic Surveillance). CorvGenSurv directly amplified viral RNA from COVID-19 patients’ Nasopharyngeal/Oropharyngeal (NP/OP) swab specimens and sequenced the SARS-CoV-2 whole genome in three segments by long-read, high-throughput sequencing. Sequencing of the whole genome in three segments significantly reduced sequencing data waste, thereby preventing dropouts in genome coverage. We validated the precision of our pipeline by both control genomic RNA sequencing and Sanger sequencing. We produced near full-length whole genome sequences from individuals who were COVID-19 test positive during April to June 2020 in Los Angeles County, California, USA. These sequences were highly diverse in the G clade with nine novel amino acid mutations including NSP12-M755I and ORF8-V117F. With its readily adaptable design, CorvGenSurv grants wide access to genomic surveillance, permitting immediate public health response to sudden threats.
Pathogen whole genome sequencing informs evidence-based public health decisions by providing crucial data for disease transmission, new outbreak detection, and vaccine candidate selection 1,2 , as demonstrated by HIV 3 , tuberculous 4 , Ebola 5,6 and Zika 7 outbreaks. In the immediate response to COVID-19, several studies have demonstrated that genomic surveillance outcomes were not only comparable to epidemiological contact tracing data 8,9,10 but also capable of tracing previously unknown linked transmissions 11 . Genomic investigation informed the public health decisions to prevent further spread of SARS-CoV-2, including travel restrictions and stay-at-home orders in response to the identification of travel-related clusters and local clusters 10 . Rapid SARS-CoV-2 whole genome sequencing is therefore an essential public health measure.
SARS-CoV-2 whole genome sequencing also has important utility in the surveillance of amino acid mutations that may result in changes to viral protein physical structures 12,13 . Such changes in structure can potentially alter virus transmissibility 14 , disease severity 15 , or reduce vaccine efficacy 16,17 . Over 90 vaccines are currently in development 18 and spike protein based vaccine candidates have induced both neutralizing antibody and T-cell responses against SARS-CoV-2 19 . However, diverse mutations in viral proteins, as observed in GISAID 20,21 and Nextstrain 22 , may lower COVID-19 vaccine effectiveness 16,17 , as observed in influenza vaccines 23,24 . Similarly, mutations in therapeutic drug target regions such as RNA-dependent RNA polymerase (RdRP) can potentially lead to treatment failure 25 . Therefore, it is important to be able to surveil virus mutations in real-time to ensure the efficacy of prophylactic vaccines and therapeutic drugs.
The greatest barriers to the deployment of massive SARS-CoV-2 whole genome sequencing for routine use are the complexity and expense of the workflows. Current SARS-CoV-2 whole genome sequencing methods, including the widely-used Artic network protocol, generally rely on short-read amplifications and sequencing 8,10,11,12,26,27,28,29,30,31 . However, a high number of short-read amplifications increases the risk of genome coverage dropouts 8 . Minimizing this risk often requires deep sequencing coverage, resulting in high-sequencing cost. To date, no SARS-CoV-2 whole genome sequencing method has concurrently accomplished high-resolution, simple workflow, and low-cost.
Herein, we introduce CorvGenSurv (Coronavirus Genomic Surveillance), an accurate and cost-efficient SARS-CoV-2 whole genome sequencing platform for COVID-19 genomic surveillance. We directly amplify viral RNA and sequence it in three segments using long-read, high-throughput sequencing. Sequencing of the SARS-CoV-2 whole genome in three segments prevents dropouts in genome coverage by minimizing ambiguous reads and eliminating assembly errors. CorvGenSurv is a streamlined, high-resolution, and cost-effective pipeline that facilitates the deployment of real-time SARS-CoV-2 whole genome sequencing for public health.
We extracted SARS-CoV-2 RNA from remnant COVID-19 positive NP/OP swab specimens in universal viral transport media. Rather than converting viral RNA to cDNA, we instead directly amplified viral RNA via long-read reverse transcriptase-polymerase chain reaction (RT-PCR). Three overlapping RT-PCR amplicons were produced for each specimen. Each specimen’s three overlapping RT-PCR amplicons were then indexed, and multiplexed specimens sequenced by long-read high-throughput sequencing 32,33 (Fig. 1a). After de-multiplexing, we built the consensus sequence of each of the three segments, correcting sequencing errors. These overlapping consensus sequences were then assembled into a near full-length SARS-CoV-2 whole genome sequence for each COVID-19 case.
We measured the accuracy of CorvGenSurv by sequencing the genomic RNA of the USA-WA1/2020 control strain (GenBank: MN985325.1) provided by BEI Resources (NR-52285). The assembled near full-length SARS-CoV-2 sequence obtained by CorvGenSurv matched 100% to the sequence of this control strain. This accuracy was achieved by our error correction step during consensus sequence building. To assess how many reads are required to produce an accurate consensus sequence, we performed a bootstrap test, randomly sampling a given number of reads and building the consensus sequence of those reads. As shown in Fig. 1b, when three reads were sampled, only 67% [52.4–78.9%] of the 1000 bootstrap runs’ resulting consensus sequences were consistent with the reference sequence of USA-WA1/2020. When we generated a consensus sequence from 31 or more reads, the consensus sequence of all 1000 bootstrap runs was identical to the control sequence. This 31-fold sequencing depth ensured the production of a SARS-CoV-2 whole genome sequence that perfectly matched the control sequence.
We also spot-checked the precision of CorvGenSurv by Sanger sequencing, the current gold standard for precision checking 34 . USA/CA-LAC-USC1 in Table 1 was selected, and a near full-length SARS-CoV-2 sequence was obtained by assembling a total of 49 overlapping segments that were sequenced by Sanger sequencing. The resulting Sanger sequence matched 100% to the corresponding USA/CA-LAC-USC1 sequence obtained by CorvGenSurv.
Figure 2b marked the amino acid changes of each of our 25 sequences in reference to the Wuhan-Hu-1 sequence. We reported a total of nine new amino acid mutations that were not observed among 28,176 global SARS-CoV-2 sequences in GISAID 20,21 . These include M755I in RdRP (NSP12) and V117F in the ORF8 protein (Table 1).
We compared the prevalence of each amino acid mutation from Wuhan-Hu-1 that had greater than 2% frequency either globally, in the USA, or in California as of July, 2020 (Fig. 2c). Mass circulation of G clade strains were confirmed by around 80% prevalence of P323L in RdRP and D614G in the S protein globally. In USA and California, the prevalence of T85I in NSP2 and Q57H in ORF3a were observed to be above 40%, as shown in Fig. 2c. The prevalence of P1263L mutation in the S protein was sixfold greater in California than the global prevalence. Additionally, S194L mutation in the N protein was much more prevalent in California, compared to other regions.
The development of a cost-effective SARS-CoV-2 genotyping protocol is crucial to expanding COVID-19 surveillance efforts. The per-specimen supply cost of our SARS-CoV-2 whole genome sequencing method was estimated to be $33.8. This includes RNA extraction from NP/OP swab media ($4.58), long-range RT-PCR ($25.8), index PCR ($2.6), and long-read high-throughput sequencing ($0.82).
We estimated the per-specimen high-throughput sequencing cost by assessing the maximum number of specimens ( \(_\) ) we can process in a single sequencing run. We observed that 31 or more reads would be required to produce an accurate consensus sequence (Fig. 1b). Therefore, we estimated how many reads are required on average to obtain a minimum of 31 reads per consensus sequence. The expected minimum value of \(_\) Poisson random variables is given by \(\sum_^<\infty><[\gamma\,\,\left(k,\lambda\right)/<\Gamma>\left(k\right)]>^<_>\) , where \(\lambda\) is the mean of the Poisson distribution, \(\gamma \left(k,\lambda \right)\) is the lower incomplete gamma function, and \(\mathrm\left(k\right)\) is the gamma function. Since we need three segments per specimen, \(_\) is given by 3 \(\times\,_\) . We obtained around 1.58 million circular consensus sequence (CCS) reads with 99% accuracy for a library size of ~ 10 kb from a single sequencing run. We estimated that 8779 specimens can be processed in a single sequencing run, where 60 reads per consensus sequence can be obtained on average to recover a minimum of 31 reads per consensus sequence. Similarly, if we aim to recover a minimum of 100 reads per consensus sequence as a conservative approach, around 3657 specimens can be processed in a single sequencing run. Assuming a sequencing cost of $3000, this yields a per specimen sequencing cost of $0.82 for large-scale sequencing efforts.
This estimated low-sequencing cost suggests a cost advantage over current short-read based SARS-CoV-2 whole genome sequencing methods 8,10,11,12,26,27,28,29,30 . The required time for CorvGenSurv is around 15 h for sample library preparation and 30 h for sequencing. Our cost-effective workflow may therefore enhance the implementation of real-time massive genomic SARS-CoV-2 surveillance.
We estimated the rate of SARS-CoV-2 evolution from the 25 genomes we obtained using a Bayesian phylogenetic reconstruction method 37,38 . The SARS-CoV-2 nucleotide evolution rate was estimated to be 7.51 \(\times\) 10 –4 substitutions per site per year [95% highest posterior density (HPD): 2.04 \(\times\) 10 –4 to 1.34 \(\times\) 10 –3 ]. The median divergence time was estimated to be December 5th, 2019 [95% HPD: December 16th, 2018–March 19th, 2020]. This estimate was 26 days prior to the first SARS-CoV-2 sequence Wuhan-Hu-1’s sample collection time, December 31st, 2019.
The rate of SARS-CoV-2 evolution was also estimated by measuring the rate of divergence from the reference Wuhan-Hu-1 sequence. We plotted our sequences’ number of nucleotide base differences from the Wuhan-Hu-1 sequence against sample collection time difference in Fig. 3. The nucleotide substitution rate in reference to Wuhan-Hu-1 was measured to be 8.62 \(\times\) 10 –4 substitutions per site per year (95% confidence interval: 7.96 \(\times\) 10 –4 to 9.24 \(\times\) 10 –4 ). This rate was consistent with the above Bayesian phylogenetic estimate.
It has been reported that several SARS-CoV-2 variants of concern can potentially escape from vaccine-induced immunity 16,17 . The emergence and circulation of such variants has prompted efforts to continuously monitor SARS-CoV-2 evolution via real-time genomic surveillance 39,40 . To further demonstrate its importance, we analyzed influenza evolution in response to vaccination. Influenza has also been shown to escape from vaccine-induced immunity, as confirmed by serologic testing 41 . The WHO changed the 2019–2020 H1N1 Northern Hemisphere Flu vaccine strain from A/Michigan/45/2015 (grey diamond in Fig. 4a–c) to A/Brisbane/02/2018 (purple diamond in Fig. 4a–c) 42 . Our maximum likelihood tree analysis showed that influenza H1NI Hemagglutinin (HA) sequences evolved away from this new vaccine strain. As shown in Fig. 4a, a total of 255 H1N1 Hemagglutinin (HA) sequences sampled in April 2019 were relatively close to the 2019–2020 H1N1 vaccine strain (purple diamond in Fig. 4a). However, we found that HA sequences collected in January 2020 (1140 sequences total) had evolved away from this vaccine strain. This rapid evolution was visualized in a two-dimensional map obtained by multidimensional scaling of the pairwise HA sequence distance matrix (Fig. 4b,c). As plotted in Fig. 4d, the nucleotide base difference from the vaccine sequence was significantly higher among the January 2020 sequences, compared to the April 2019 sequences (median distance 18 vs. 23, p < 0.001).
This observed evolution away from the vaccine strain may explain the 2019–2020 seasonal influenza vaccine’s low effectiveness against influenza A (H1N1). Influenza vaccine efficacy was estimated to be 37% against H1N1 strains during the 2019–20 flu season 43 . Viral evolution in response to vaccination remains a challenge for vaccine design, and the 2019–20 flu season highlights the need for real-time surveillance of viral sequences. Like influenza, SARS-CoV-2 has the potential to evolve in response to vaccination 44 , and thus SARS-CoV-2 evolution must be continuously monitored in real time to prevent vaccine failure and guide future vaccine strain selection.
We developed CorvGenSurv (Coronavirus Genomic Surveillance), a streamlined, high-resolution, and cost-effective SARS-CoV-2 whole genome sequencing method in which viral RNA is directly amplified, indexed, and sequenced. CorvGenSurv’s precision and accuracy have been validated by Sanger sequencing and control genomic RNA sequencing, respectively. Our targeted amplification of the whole genome in three segments minimized the risk of genome coverage dropouts. Additionally, our streamlined long-read sequencing protocol significantly reduced workflow complexity and thus has a potential cost advantage over currently available short-read sequencing approaches 8,10,11,12,26,27,28,31 .
We obtained 25 whole genome sequences from NP/OP remnant COVID-19 positive test specimens that were collected from April to June 2020 at Los Angeles County, California, USA. The Los Angeles County pandemic accounted for 33% of all COVID-19 cases in California with around 3500 new cases per day in July 2020 45,46 . Our maximum likelihood tree analysis showed that the 25 strains were highly diverse within the G, GR and GH clades 35 . The G clade and tis lineage have been dominant since late March, 2020 and its D614G mutation in the spike protein has been associated with increased transmissibility 14 . In addition to the G clade’s other common mutation, NSP12-P323L, frequent mutations we observed include N-S194L, N-R203K, N-G204R, NSP2-T85I, and ORF3a-Q57H. From this wide spectrum of viral mutations, we estimated SARS-CoV-2 evolutionary rate by a Bayesian phylogenetic reconstruction method 37,38 and divergence measures. These two estimates were 7.51 \(\times\) 10 –4 and 8.62 \(\times\) 10 –4 substitutions per site per year respectively, which were consistent with recent reports on the SARS-CoV-2 evolutionary rate 47,48,49 .
We identified new amino acid mutations in the NSP3, NSP4, RdRP (NSP12), ORF8, and N proteins. Mutations in these proteins have the potential to trigger immune escape, alter viral replication capacity, or modulate viral immune suppression. Numerous T cell epitopes have been annotated in these proteins and thus mutations can lead to viral escape from cellular immune responses 50,51 . While mutations in the N, NSP3, NSP4 or NSP12 proteins can alter viral replication capacity 13,52,53,54,55 , those in the ORF8 or N proteins may dysregulate host innate immune responses via type I interferon signal modulation 54 . It is therefore of great importance to monitor mutations from a diverse array of proteins by whole genome sequencing to pinpoint those that are relevant to viral pathogenicity and immune escape.
SARS-CoV-2’s S protein is the key target for many of the vaccines currently in development 18 and thus mutations in this region are being closely monitored for their potential to render these vaccines ineffective. The observed influenza HA sequence evolution away from this year’s vaccine strain clearly indicates the necessity of screening mutations in viral surface proteins. Genomic surveillance data can also inform SARS-CoV-2 drug resistance. Mutations in RdRP (NSP12) can potentially result in decreased binding affinity with nucleoside analogs such as Remdesivir and Favipiravir 56 . Therefore, genomic surveillance on circulating pandemic strains is crucial for guiding strategies for COVID-19 vaccination and therapeutic intervention and our simple and cost-efficient pipeline has the potential to advance real-time COVID-19 surveillance.
CorvGenSurv can accurately survey new mutations in real-time and thus provide crucial data for disease transmission, new outbreak detection, and vaccine candidate selection. As an alternative to the more common short-read sequencing approaches 8,10,11,12,26,27,28,29,30 , CorvGenSurv can facilitate the widespread adoption of real-time COVID-19 genomic surveillance and can therefore serve as an important tool for public health decision-making.
We accessed 25 de-identified nasopharyngeal (NP) and oropharyngeal (OP) remnant specimens which tested positive for COVID-19 via the Roche cobas ® SARS-CoV-2 qualitative EUA assay at USC Clinical Laboratories, Keck Medicine of USC. These specimens were collected from Los Angeles County, California, USA between April 13th and June 22nd, 2020 (Table 1). This study (HS-20-00326) was approved by the Institutional Review Board of the University of Southern California as non-human subjects research. Table 1 lists each specimen’s collection date and cycle threshold (Ct) values for the ORF1 a/b non-structural SARS-CoV-2 unique region (Ct-1) and the pan-Sarbecovirus conserved region in the structural protein envelope E-gene (Ct-2).
We downloaded near full-length SARS-CoV-2 sequences from GISAID that were registered by July 27th, 2020. The sequences were globally aligned using MUSCLE 57 , and their 5’- and 3’-ends trimmed-out, yielding a 597–29340 segment of the Wuhan-Hu-1 reference strain. We removed any sequences that were shorter than this region and sequences with one or more ambiguous bases to obtain a total of 28,176 near full-length SARS-CoV-2 sequences, including 9339 originating from the USA and 1215 from California. In Table 1, all amino acid mutations from the reference sequence were annotated and novel mutations in our 25 sequences were reported.
SARS-CoV-2 RNA was extracted from remnant NP/OP swab specimens in universal transport media using a MagMAX™ Viral/Pathogen Nucleic Acid Isolation Kit with a KingFisher Duo Prime automated nucleic acid purification system (Thermo Fisher Scientific). A total of 20 μl (15 \(\le\) Ct \(t \(t ≥ 25) of viral media was diluted in 1× PBS buffer to a total volume of 400 μl and added to a deep-well plate for extraction using the manufacturer’s protocol. Around 80 μl of extracted SARS-CoV-2 RNA was recovered. This RNA was used as input for three separate 10 kb long RT-PCRs using the SuperScript™ IV One-Step RT-PCR System (Thermo Fisher Scientific). As listed in Supplementary Table S1, the RT-PCR primers were 1_LEFT/33_RIGHT or For-c/Rev-c for the first segment, For-d/67_RIGHT or 33_LEFT/67_RIGHT for the second segment, and 67_LEFT/Rev-e or 67_LEFT/98_RIGHT for the third segment. Here, 1_LEFT, 33_RIGHT, 33_LEFT, 67_RIGHT, 67_LEFT, and 98_RIGHT were the ARTIC network primers (V3, https://artic.network/ncov-2019), where the "nCoV-2019_" suffix was removed for simplicity. Primers, For-c, Rev-c, For-d, and Rev-e were designed in-house. These primers span the near full-length SARS-CoV-2 genome in three segments with overlap for genome assembly. The master mix was composed of 25 μl RT-PCR master mix, 2.5 μl each of 10 μM forward and reverse primers, 0.5 μl SuperScript IV RT mix, and 19.5 μl of viral RNA to a total volume of 50 μl. The samples were cycled according to the manufacturer’s instructions: 50 °C for 10 min of RT, followed by 98 °C for 2 min for RT inactivation, then 35 cycles of 98 °C for 10 s, 55 °C for 10 s, 72 °C for 5 min, then a final extension of 72 °C for 5 min, held at 4 °C until the next step.
The RT-PCR products were then purified and concentrated with Ampure XP beads (Beckman Coulter) with a 0.8 × bead volume, two 70% ethanol washes, and elution in 25 μl Nuclease-Free H2O. These purified RT-PCR products were subjected to long-range index PCR using PrimeSTAR GXL DNA Polymerase (Takara Bio) according to the manufacturer’s instructions. After index PCR, product bands were confirmed via 1% agarose gel electrophoresis with the E-Gel Electrophoresis System (ThermoFisher, CA). Samples were equimolar pooled and shipped to DNA Technologies and Expression Analysis Core at UC Davis Genome Center for long-read, single-molecule, real-time (SMRT) sequencing on the PacBio Sequel II system with a 30-h movie.
Specimen, USA/CA-LAC-USC1, was selected for near full-length SARS-CoV-2 sequence confirmation via Sanger sequencing. Purified RT-PCR product for each of the three segments of USA/CA-LAC-USC1 was amplified via long-range PCR with the master mix composed of 10 μl of PrimeSTAR GXL buffer, 4 μl of dNTPs, 2.5 μl each of 2 μM forward and reverse primers, 1 μl of PrimeSTAR GXL DNA polymerase, 26 μl nuclease free H2O, and 4.0 μl of purified RT-PCR product to a total volume of 50 μl. The samples PCR cycled with 25 cycles of 98 °C for 10 s and 68 °C for 10 min, held at 4 °C until the next step. These products were Ampure cleaned with a 0.8 × bead volume and diluted to 20 ng/μl. A total of 10 μl of the diluted product was combined with 5 μl of each Sanger sequencing primer in 5 μM. As listed in Supplementary Table S2, a total of 49 primers were selected from the V3 ARTIC network primers (https://artic.network/ncov-2019) with minor modifications to minimize self-dimer and hairpin formation. The samples were then shipped to GENEWIZ (South Plainfield, NJ) for Sanger sequencing. The obtained sequences were subject to contig assembly using the CAP3 Sequence Assembly 58 followed by manual inspection, yielding USA/CA-LAC-USC1’s near full-length SARS-CoV-2 genome (55–29,836).
Each SARS-CoV-2 genome segment was de-multiplexed based on their indexes, creating one fasta file for each of the three segments per specimen. 100 reads were randomly picked from each fasta file and aligned using MUSCLE 57 to obtain a consensus sequence for each segment. The three segment consensus sequences were then assembled to produce each specimen’s near full-length SARS-CoV-2 whole genome sequence (405–29,431).
We used the PHYML program to produce maximum likelihood trees 59 . The general time-reversible model was used with the ‘ML’ option, and invariable sites were estimated with 12 substitution rate categories. The tree was generated by the BIONJ option. We used FigTree to present the obtained tree file.
BEAST42 v.1.10.4. was used to estimate the rate of SARS-CoV-2 evolution and divergence time 37 . The 25 SARS-CoV-2 near full-length genome sequences obtained in this study were input to BEAST along with its sample collection time (Table 1). Uncorrelated relaxed clock model with log-normal distribution was assumed with flexible skygrid coalescent tree priors and a single GTR + Γ substitution model 38 . The Monte Carlo (MCMC) length of chain was 2 \(\times\) 10 6 . The median evolution rate and divergence time along with respective 95% highest posterior density (HPD) were reported.
We globally aligned 5772 complete HA genome sequences that were collected between April 2019 and March 2020 along with two H1N1 vaccine strains A/Michigan/45/2015 and A/Brisbane/02/2018. These sequences were downloaded from GISIAD 20,21 . The pairwise nucleotide base differences, Hamming Distance (HD), was measured and this HD matrix was input to a multi-dimensional scaling (MDS) algorithm from the python scikit-learn package (0.19.1). This MDS algorithm converted the HD matrix into a two-dimensional representation of the HA sequences by minimizing stress, the sum of squared distance of the disparities. Each sequence’s position was updated until the stress difference was less than 0.1, and by repeating this procedure 400 times, the configuration with the minimum stress value was selected.
SARS-CoV-2 whole genome sequences sequenced in this study are available in GISAID (accession numbers EPI_ISL_569664-EIP_ISL_569688).
We thank the Molecular Pathology Division of the USC Keck Clinical Laboratories for providing the remnant NP/OP specimens that are sequenced in this study. The genomic RNA of USA-WA1/2020 strain was deposited by the Centers for Disease Control and Prevention and obtained through BEI Resources, NIAID, NIH: Genomic RNA from SARS-Related Coronavirus 2, Isolate USA-WA1/2020, NR-52285. We thank Dr. Jing-Hsiung James Ou, Dr. Sarah Hamm-Alvarez, and Dr. Thomas Buchanan for their help. We thank all laboratories and authors who shared global SARS-CoV-2 whole genome sequences and influenza H1NI Hemagglutinin (HA) sequences by depositing to GISAID’s EpiFlu (TM) Database.
