key: cord-0913625-rslt05hq
authors: Rosenthal, Sun Hee; Gerasimova, Anna; Ruiz-Vega, Rolando; Livingston, Kayla; Kagan, Ron M.; Liu, Yan; Anderson, Ben; Owen, Renius; Bernstein, Laurence; Smolgovsky, Alla; Xu, Dong; Chen, Rebecca; Grupe, Andrew; Tanpaiboon, Pranoot; Lacbawan, Felicitas
title: Development and validation of a high throughput SARS-CoV-2 whole genome sequencing workflow in a clinical laboratory
date: 2022-02-08
journal: Sci Rep
DOI: 10.1038/s41598-022-06091-0
sha: 4a2e611e80a93b283a677d2287b5f5061039cf1c
doc_id: 913625
cord_uid: rslt05hq

Monitoring new mutations in SARS-CoV-2 provides crucial information for identifying diagnostic and therapeutic targets and important insights to achieve a more effective COVID-19 control strategy. Next generation sequencing (NGS) technologies have been widely used for whole genome sequencing (WGS) of SARS-CoV-2. While various NGS methods have been reported, one chief limitation has been the complexity of the workflow, limiting the scalability. Here, we overcome this limitation by designing a laboratory workflow optimized for high-throughput studies. The workflow utilizes modified ARTIC network v3 primers for SARS-CoV-2 whole genome amplification. NGS libraries were prepared by a 2-step PCR method, similar to a previously reported tailed PCR method, with further optimizations to improve amplicon balance, to minimize amplicon dropout for viral genomes harboring primer-binding site mutation(s), and to integrate robotic liquid handlers. Validation studies demonstrated that the optimized workflow can process up to 2688 samples in a single sequencing run without compromising sensitivity and accuracy and with fewer amplicon dropout events compared to the standard ARTIC protocol. We additionally report results for over 65,000 SARS-CoV-2 whole genome sequences from clinical specimens collected in the United States between January and September of 2021, as part of an ongoing national genomics surveillance effort.

High-throughput NGS workflow optimization. Our laboratory has developed an automated, highthroughput SARS-CoV-2 NGS workflow (Fig. 1) . This workflow utilizes a 2-step PCR NGS library preparation method: (1) gene-specific PCR to amplify the SARS-CoV-2 whole genome, using primers published by the ARTIC network, with modifications to add Illumina sequencing primer binding sites; and (2) index PCR to add specimen-specific barcoded sequencing adapters by fusion PCR using the Illumina sequencing primer binding sites. We first optimized primer pools to give even coverage across the SARS-CoV-2 genome. Gene specific primers were pooled into 4 pools (pool 1A, 1B, 2A, 2B, Supplementary Table S1), and positive patient specimens with RT-qPCR cycle threshold (Ct) value of 24 (n = 11) were tested. To evaluate the coverage uniformity, we computed the average coverage of each amplicon at a normalized depth of 200,000 mapped reads. All 98 amplicons produced adequate coverage with a mean coverage of 973X (SD, 719; CV, 73.9%) ( Fig. 2A) . However, 7 amplicons, including 3 in the spike protein coding region, produced relatively lower coverage. When the same sample www.nature.com/scientificreports/ set was tested with a standard ARTIC v3 method, even coverage across the whole genome was achieved with a mean coverage of 1390X (SD, 658; CV, 47.3%) at 200,000 mapped reads (Fig. 2B ). Next, to improve coverage for these low-performing amplicons, we increased the corresponding primer concentration in the same primer pool or added the primer set to a different pool. When tested with RT-qPCR positive patient specimens (Ct 24, n = 5), the optimized primer pools improved the coverage of these low-performing amplicons by 2-to fivefold. The coverage for amplicons 26, 68, 74 and 75 was improved approximately twofold from 106X, 52X, 61X and 130X to 191X, 273X, 100X and 316X, respectively. In addition, the coverage for amplicons 29, 67 and 79 was improved approximately 4-to fivefold from 165X, 52X and 190X to 773X, 186X, and 720X, respectively. These changes resulted in improved amplicon balance with a mean coverage of 893X (SD, 514; CV, 57.6%) (Fig. 2C) . For the same sample set, the standard ARTIC v3 method achieved 1470X mean coverage depth (SD, 692; CV 47.1%) (Fig. 2D) .

As the SARS-CoV-2 genome has evolved, various mutations were found at the ARTIC v3 primer binding sites. Upon sequence analysis of 3506 positive samples collected and processed in January and February of 2021, we found that 96% (3367/3506) of the samples had at least one primer binding site mutation with a median of 3 (range: 1-9) mutations with the potential to impair PCR efficiency per affected sample (Supplementary Table S3) .

To minimize adverse effects on sequencing coverage, we employed a modified touchdown PCR method by gradually reducing the annealing temperature from 65 °C to 55 °C (0.7 °C/s) within each PCR cycle. We compared 79 samples with two different annealing temperature settings. With this approach, the percent amplicon dropout due to primer binding site mutations was decreased from 0.50% to 0.01% (Table 1 ). The percent amplicon dropout was calculated from the number of amplicons that did not generate coverage divided by the total amplicon number. Each dropout amplicon was manually reviewed for the presence of a mutation at the affected primer binding site.

We assessed intra-assay precision using 188 unique specimens with Ct values between 10 and 25 run in 3 replicates. Most (523/564; 92.8%) sequences met quality control metrics (Fig. 3A) ; 175/188 (93.1%) unique samples met quality control requirements for 2 or more replicates and were utilized for lineage and clade comparison. All 175 samples had concordant clade and lineage assignments for all replicates (Table 2) . Fig. S1 A) . The average %CV for variant frequency between replicates was 0.5% (min CV: 0.0%; max CV: 6.7%). www.nature.com/scientificreports/ We then assessed inter-assay precision for 168 unique specimens with Ct values between 10 and 25 over 3 independent runs. Overall, 479/504 (95.0%) sequences met quality control requirements (Fig. 3B ). Over 95% (160/168) of samples had valid sequence data for 2 or more replicates and were utilized for clade and lineage comparisons. All 160 samples had concordant clade and lineage assignments ( Table 2 ). On average, there were 5.9 ± 3.2 spike protein amino acid substitutions per sample with 99.5% (946/951, 95% CI: 98.8-99.8%) qualitative agreement between runs when the minority variant frequency was present in > 20% of the reads ( Supplementary  Fig. S1 B) . The average %CV for the variant frequency between replicates was 0.9% (min CV: 0.0%; max CV: 20.3%).

Assay sensitivity. To demonstrate that our high-throughput workflow offers adequate sensitivity to yield complete SARS-CoV-2 genomes, we determined the limit of detection. The limit of detection was defined as Table 2 . Intra-and inter-assay concordance of clade and lineage assignment using the automated, highthroughput SARS-CoV-2 NGS workflow. a SARS-CoV-2 clades were assigned with Nextclade v1.3.0 (https:// clades. nexts train. org/) and lineages were assigned with Pangolin v3.1.11 (https:// pango lin. cog-uk. io/). Accuracy study. Next, we assessed the accuracy of the consensus sequences generated by this workflow.

Three commercial synthetic RNA positive controls (clade Alpha, Beta, and Gamma variants; Twist Bioscience) were tested 12 times each, using 4 replicates per set-up in 3 independent set-ups, by 3 different scientists. Of note, the synthetic RNA controls do not cover 100% of the SARS-CoV-2 genome, and amplicon dropout was expected. In 36 trials, the mean percent consensus sequence was 95.1% (min 93%; max 98%). All controls resulted in 100% positive percent agreement for clade and lineage assignments (95% CI: 90.4%, 100.0%) and 100% for spike protein mutations that were detected when the minor allele frequency was over 20% (95% CI: 98.8%, 100.0%).

In addition, a total of 84 qRT-PCR negative samples were evaluated, in 3 independent set-ups using 28 negative samples per run, by 3 different scientists. All negative samples gave negative results, with mean amplicon coverage of 0.2 (min 0.0%, max 0.6%) relative to the PCR plate average, yielding 100% sample-level negative percent agreement (95% CI: 95.6%, 100.0%).

To assess assay robustness, a total of 2688 samples, including 1662 unique SARS-CoV-2 positive clinical specimens, 21 synthetic positive controls, 49 SARS-CoV-2 negative specimens and 56 no-template controls (NTCs), were sequenced in a single Illumina NovaSeq sequencing run. Over 94.3% (2416/2562) of the total number of sequenced samples and 95.3% (1584/1662) of the unique clinical specimens passed quality control with a median coverage of 2478X. All negative specimens gave a negative result (median coverage Tables S4, S5 ). Moreover, there was 100% clade and lineage concordance for all positive clinical specimens processed in the same 384-well plates with the 2 higher-coverage NTC wells. These results indicate that a low degree of NTC contamination is unlikely to interfere with accurate consensus sequence generation or with clade and lineage assignments. Next, we analyzed the proportion of samples that generated 100% consensus coverage out of the samples that passed coverage QC (Table 3) . With the optimized high-throughput workflow, over 95.9% of samples generated a complete consensus (range: 83.0% to 100% per clade). When the same set of samples were analyzed with a standard ARTIC v3 workflow, the proportion of samples with complete consensus sequence was much lower, 66.1% (range 0% to 100.0% per clade). In some variants sequenced with the ARTIC v3 workflow (clades Beta, Delta, Epsilon, and Lambda), amplicon dropouts were observed in almost all cases; dropout was resolved by employing our high-throughput workflow utilizing touchdown PCR.

Large-scale surveillance study. From January to September 2021, over 65,000 SARS-CoV-2 qRT-PCRpositive or transcription-mediated amplification (TMA)-positive specimens were successfully sequenced in our laboratory on behalf of the CDC National SARS-CoV-2 Strain Surveillance (NS3) project (https:// www. aphl. org/ progr ams/ prepa redne ss/ Crisis-Manag ement/ COVID-19-Respo nse/ Pages/ Seque nce-Based-Surve illan ce-Submi ssion. aspx). Specimens were collected from all 50 US states, the District of Columbia and Puerto Rico. The leading 10 states that cumulatively contributed two-thirds (66%) of the specimens over this study period were Texas (15.1%), California (11.8%), Florida (10.5%), Pennsylvania (6.1%), Nevada (5.8%), Massachusetts (4.4%), Illinois (4.1%), Connecticut (4.1%) and Ohio (4.0%). Between January and the end of May 2021, the sequencing libraries were constructed using the modified ARTIC v3 protocol as described in "Methods". Subsequently, libraries were constructed using the high-throughput workflow designed to improve scalability, as reported in this study. Cumulatively, the Delta variant of concern (VOC) was detected at the highest frequency (35%), because of the steep increase in prevalence of this variant after June 2021 (https:// covid. cdc. gov/ covid-datatrack er/# varia nt-propo rtions). Viruses in clade Alpha were the second-most prevalent, accounting for 29% of specimens tested, reflecting their high prevalence before June 2021. Other variants being monitored (VBMs) were detected at much lower frequencies, including Epsilon (5%), Iota (5%), and Gamma (4%) ( Supplementary  Table S8 ). Beta, Lambda, Eta, and Kappa variants were found in less than 1% of the specimens sequenced.

We tracked the proportions of VBM and VOC variants on a weekly basis (Fig. 5) . The Epsilon variant was predominant (25.6%) in January and February 2021, then subsequently declined in frequency as the Alpha and Iota variants became more prevalent (65.3% and 13.3%, respectively) through May 2021. The Delta variant 

Pathogen WGS has been employed to investigate cases of foodborne disease outbreaks, methicillin-resistant Staphylococcus aureus (MRSA) outbreaks, and outbreaks of other bacterial and fungal pathogens 20 . However, these studies were retrospective in nature. Rapid high-throughput viral genomic sequencing is being increasingly used for outbreak research and is changing the practice of epidemiology 21 . For example, advances in high-throughput sequencing facilitated near-real time WGS to investigate the 2014-2015 Ebola virus outbreak in West Africa 22 . In the present COVID-19 pandemic, rapid SARS-CoV-2 WGS, combined with epidemiologic methods, identified or excluded nosocomial transmission events to directly affect outbreak management in real time 23 . Hospital-based WGS set up in collaboration with public health departments has also facilitated real-time outbreak investigations and identification of transmission chains 24 . At a national level, the Coronavirus Disease 2019 (COVID-19) Genomics UK Consortium (COG-UK) has developed high-throughput sequencing and analysis workflows to generate hundreds of thousands of SARS-CoV-2 genomic sequences, enabling large-scale genomic epidemiology to inform public health response (https:// www. cogco nsort ium. uk/). In the United States, the CDC established the SARS-CoV-2 Sequencing for Public Health Emergency Response, Epidemiology and Surveillance (SPHERES) initiative in the early stages of the pandemic to accelerate the use of real-time pathogen sequence data and molecular epidemiology for the COVID-19 pandemic response (https:// www. cdc. gov/ coron avirus/ 2019-ncov/ varia nts/ spher es. html). Large scale and high-throughput WGS of SARS-CoV-2 is essential to enable real-time epidemiologic surveillance and to better understand the evolution and spread of the pandemic, improve epidemiologic tracking, and enhance treatment and prevention strategies for COVID-19. To date, several groups have developed different protocols and workflows for SARS-CoV-2 WGS by NGS. However, most of the published studies are small in study sample size, low in throughput or may require complex enzymatic steps 11, 12, [25] [26] [27] [28] [29] . In one study, the commercially available COVIDSeq test (Illumina) was scaled up to sequence 752 clinical samples in duplicate for a total of 1536 samples and controls on a NovaSeq 6000 instrument using an S4 flow cell 30 . However, this study described only a single sequencing run and did not utilize an automated high throughput workflow. In the current study, we were able to use the lower output and less expensive SP flow cell to sequence over 2600 samples.

Here, we report an automated, high-throughput workflow for SARS-CoV-2 WGS optimized for large scale surveillance studies, utilizing a 2-step PCR NGS library preparation method with modified ARTIC protocol v3 primers. Validation studies performed on 1711 unique clinical samples demonstrated high precision (100% interand intra-assay precision) and accuracy (100% PPA and 100% NPA). Slightly reduced, but adequate sensitivity was achieved in comparison to the lower-throughput in-house validated ARTIC v3 protocol: near complete consensus sequence was generated for samples of qRT-PCR Ct values less than 28 using the 2-step PCR workflow (Fig. 4) and with samples of Ct less than 30 by the ARTIC v3 workflow (data not shown). www.nature.com/scientificreports/ For large-scale surveillance studies, integration of robotic liquid handlers is a crucial component. Various automation platforms are available for NGS library preparation. Upon evaluation of frequently utilized liquid handlers, the Agilent Bravo platform produced the most consistent results for our workflow in a 384-well format and was used for library preparation. As noted by Gohl et al. 17 , a 2-step PCR method increased primer dimer formation which can result in lower sequencing quality. To remove any primer dimer present, we implemented 3 library cleaning steps: (1) Ampure bead clean up using a centrifugal washer in 384-well format; (2) manual Ampure bead clean up of a pooled library; (3) size selection using 1.5% gel on a BluePippin station prior to loading on the Illumina NovaSeq 6000. With these stringent clean-up steps, our validation runs achieved over 90% (range 90.6%-93.4%) high base call Phred quality scores (Q30 or higher).

The ARTIC amplicon-based whole-genome amplification of SARS-CoV-2 is considered a highly sensitive and accurate method and is currently employed by many laboratories around the world. One problem associated with amplification of the SARS-CoV-2 genome is amplicon dropout which can be caused when a given primer cannot stably hybridize to its specific complementary sequence binding site because of novel mutations at the primer binding site. As the SARS-CoV-2 genome has evolved, various mutations were found in ARTIC v3 primer binding sequence 31 . This type of dropout often requires redesign and revalidation of primers. To minimize amplicon dropout, we employed a modified touchdown PCR method which reduced amplicon dropout ( Table 1 ). The touchdown PCR method is expected to increase PCR specificity by reducing the amplification of spurious amplicons at higher annealing temperatures, while at the same time maintain a competitive advantage for the targeted amplicons and increase their yield, even in the presence of potential primer-binding site mutations, as the annealing temperature is gradually reduced 19 . The effects of this improvement were most prominent for Beta, Delta, Epsilon, and Lambda clade specimens, achieving 100% consensus sequence coverage in most trials (Table 3) . On the other hand, for these 4 clades, when analyzed by the standard ARTIC v3 workflow, nearly all samples failed to achieve full consensus sequence coverage due to primer binding site mutations causing amplicon dropout (Table 3, Supplementary Table S6) . We extended our analysis to samples collected between May and September of 2021 to show that the benefit of using the touchdown PCR method is not limited to the validation sample set. Similar results were observed for Beta, Delta, Epsilon, and Lambda clade specimens with 95% of samples generating a complete consensus sequence using the high-throughput workflow whereas only 2.9% of samples generated complete consensus sequence by the ARTIC v3 workflow (Supplementary Table S7 ). The ARTIC network has published a newer ARTIC v4 primer set, which was designed to avoid current highfrequency mutations, with the goals of minimizing amplicon dropout and producing high-accuracy variant calling. However, analytical validation is required prior to using this new primer set in surveillance studies. We envision that a touchdown PCR method can be easily adapted to any workflow and can avoid frequent primer redesign and validation.

The 2-step PCR method generated 30% lower coverage compared to the ARTIC v3 method when the per sample coverage was averaged across all amplicons and normalized to 200,000 mapped reads per sample ( Fig. 2A,B) . We postulate that this difference may be due to the higher primer dimer formation during the 2-step PCR as mentioned by Gohl et al. 17 . Although the normalized coverage was reduced by approximately 30%, on average, for the 2-step PCR method, the normalized coverage still exceeded the minimum coverage necessary to reliably generate a viral genome consensus sequence or to detect viral variants present in at least 50% of the reads (consensus calling threshold). Indeed, the absolute coverage obtained even when processing > 2600 samples on the NovaSeq sequencer (median coverage 2478X) greatly exceeded the absolute coverage obtained for ARTIC v3 runs on the MiSeq sequencer (median coverage 1098X with a batch size of 192). Therefore, these differences in normalized coverage have minimal practical impact on the ability of the 2-step PCR method to generate reliable whole genome consensus sequences or to detect variants.

In summary, building on SARS-CoV-2 WGS methodological improvements reported elsewhere 17 , we report on the development and the validation of an automated 2-step PCR, high-throughput NGS workflow for SARS-CoV-2 WGS with high sensitivity and accuracy with much less amplicon dropout, which we have implemented for high-throughput SARS-CoV-2 genomic surveillance. The combination of automation and optimized bench and analysis workflows has enabled efficient, large-scale SARS-CoV-2 surveillance studies.

Sample collection. We obtained remnant RNA stored at − 80 °C, extracted from deidentified clinical specimens previously tested for SARS-CoV-2 by qRT-PCR at Quest Diagnostics between February and August 2021. For workflow development, 16 positive samples with a Cycle threshold (Ct) value 24 were processed for sequencing. Analytical method validation studies were performed on a random set of 1711 unique clinical specimens collected in March and April of 2021. In addition, 3 twist Synthetic SARS-CoV-2 RNA controls (Twist Bioscience, San Francisco, CA) and 24 negative samples (Ct > 50 by qRT-PCR) were included for accuracy studies.

For surveillance studies, a random set of residual nasopharyngeal or oropharyngeal swab specimens from SARS-CoV-2 qRT-PCR positive or transcription mediated amplification (TMA) positive specimens, performed at Quest Diagnostics, were collected across the United States between January and September 2021. For SARS-CoV-2 qRT-PCR positive sample selection, Ct < 25 was applied for qualification. In accordance with ethical requirements and in accordance with US Department of Health and Human Services guidelines for the use of deidentified specimen remnants in research studies, specimens were deidentified prior to the study and were limited to remnant extracted RNA from discarded specimens previously submitted for commercial testing. Table S1 ). The primers were based on the ARTIC v3 primer set 32 with modifications to add Illumina sequencing primer binding sites: 5'-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT-3' was added to 5' of all forward primers; 5'-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATCT-3' was added to 5' of all reverse primers. To avoid primer interactions, 11 of the ARTIC v3 primers were replaced with primers reported by Itokawa et al. 33 : forward primers for amplicons 21, 36, and 89; reverse primers for amplicons 1, 9, 13, 15, 32, 38, 76 ARTIC v3 library preparation. The cDNA synthesis and PCR steps were identical to the original ARTIC v3 protocol (https:// www. proto cols. io/ view/ ncov-2019-seque ncing-proto col-bbmui k6w) with minor modifications. In brief, RNA samples were reverse transcribed using the SuperScript IV kit (ThermoFisher, Waltham, MA) and random hexamers. RNA samples with Ct < 18 were diluted 50-fold (Ct 12-15) or tenfold (Ct 15-18) before use. In brief, 11 μl RNA samples were mixed with 1 μl 50 μM random hexamers and 1 μl 10 mM dNTPs. The mixture was incubated for 5 min at 65 °C and placed directly on ice. After 2-3 min, 7 μl enzyme mix containing 4 μl 5X Buffer, 1 μl 0.1 M DTT, 1 μl RNaseOUT RNase inhibitor, and 1 μl SuperScript IV reverse transcriptase was added to the samples. The reactions were incubated at 42 °C for 50 min, and at 70 °C for 10 min followed by cooling to 4 °C. For the multiplex PCR reactions, 2.5 μl of the cDNA was used in 25 μl reaction mixture of Q5 Hot START DNA Polymerase kit (New England Biolabs, Ipswich, MA): 5 μl 5X buffer, 0.5 μl 10 mM dNTPs, 0.25 μl polymerase and 3.6 μl 10 μM primer pool 1 or pool 2, and 13.15 μl nuclease-free water. The primer pools were based on the ARTIC v3 primer pool scheme with further optimization to enhance low-performing amplicons (Supplementary Table S2 ). The thermal program was identical to the original ARTIC protocol: 30 s polymerase activation at 98 °C followed by 30 cycles of 15 s denaturing at 98 °C and 5 min annealing and extension at 65 °C. The PCR products in pool 1 and 2 reactions for same clinical samples were combined and purified with 1X AMPure XP and quantified with Qubit Broad Range Kit on SpectraMax (Molecular Devices, San Jose, CA).

The purified PCR products were converted to Illumina sequencer-compatible libraries using the Twist Library Preparation kit (Twist Bioscience, San Francisco, CA). In brief, 300 ng of purified PCR products in 15 μl were end-repaired in 25 μl of end repair reaction mixture: 2.5 μl 10X buffer, 5 μl enzyme mix, and 2.5 μl water and incubated at 20 °C for 30 min, and at 65 °C for 30 min. Illumina adapter was ligated to the end-repaired mix by adding 25 μl ligation master mix: 10 μl 5X buffer, 5 μl enzyme mix, 2.5 μl 5 μM adapter, and 7.5 μl water and incubated at 20 °C for 15 min. The ligated products were purified using 0.8X Ampure XP beads for subsequent index PCR using KAPA HiFi HotStart ReadyMix in 50 μl reaction volume with 6 μl 5 μM index primers using the following program: 98 °C for 45 s followed by 8 cycles of 98 °C for 15 s, 60 °C for 30 s, and 72 °C for 30 s, and the final extension at 72 °C for 1 min. The indexed products were purified using 1X Ampure XP beads, quantified using Qubit Broad range reagent, and normalized by pooling 400 ng per sample. The pooled library was used for MiSeq loading with 192 samples per run.

Sequencing. For sequencing on a NovaSeq 6000 (Illumina, San Diego, CA), the final library was diluted to 1700 pM and PhiX was spiked at 18%. The combined library was denatured with 0.2 N NaOH and neutralized with 0.4 N Tris-HCl and sequenced on a NovaSeq 6000 SP Reagent kit using 2 × 251 cycles. For sequencing on a MiSeq (Illumina, San Diego, CA), the final library was diluted to 2 nM, and spiked in 20% PhiX, denatured with 0.2 N NaOH and neutralized with 0.2 N HCl, and diluted with Illumina HT1 buffer to 10 pM and sequenced using a MiSeq 600 cycle v3 kit, 2 × 251 cycles. www.nature.com/scientificreports/ were down-sampled by 50%. The fastq files were mapped to SARS-CoV-2 (MN908947.3) supplemented with the human genome (GRCh37) reference sequence using BWA v0.7.17-r1188 ( Supplementary Fig. S2 ) 34 . PCR primer sequences were trimmed from the mapped reads to MN908947.3 using iVar v1.3.1 or version 1.2.2 prior to February 2021 ( Supplementary Fig. S2 ) 35 . iVar was also used for variant calling and consensus sequence creation with minimum sequencing depth set to 10 reads for MiSeq data, and a minimum variant frequency of 0.5. For NovaSeq data, coverage equivalent to the number of mapped reads divided by 20,000 was required for consensus sequence generation and variant calling. Reads sorting and filtering were performed by SAMtools v1.9 ( Supplementary Fig. S2 ) 36 . Coverage depth of each nucleotide position was determined with BEDTools v2.29.2 ( Supplementary Fig. S2 ) 37 . Percent genome coverage was determined by counting the number of nucleotides meeting the minimum coverage requirement divided by the total SARS-CoV-2 genome length excluding the 5' and 3' ends not covered by the amplicon panel. For clade and lineage assignment, Nextclade version 1.3.0 and Pangolin version 3.1.11 with pangoLEARN 2021-08-24 were used. Clade and lineage assignments initially performed with earlier versions of Nextclade and Pangolin were subsequently reprocessed with the above-indicated versions.

Analytical method validation. Using the optimized workflow, we performed analytical method validation studies on 1711 unique clinical specimens collected in March and April of 2021, along with positive and negative controls. A sequencing batch size of 1536 or 2688 was employed using 384-well plates (4 plates for the 1536 batch size and 7 plates for the 2688 batch size). When the batch size was 1536, the resulting fastq files were down-sampled by 50% before analysis. For sample quality control, we set a minimum sequencing coverage threshold of 10X spanning at least 97% SARS-CoV-2 genome.

A novel coronavirus from patients with pneumonia in China

A new coronavirus associated with human respiratory disease in China

Novel SARS-CoV-2 variants: the pandemics within the pandemic

SARS-CoV-2 variants, spike mutations and immune escape

The biological and clinical significance of emerging SARS-CoV-2 variants

Unraveling the zoonotic origin and transmission of SARS-CoV-2

An implementation science approach to evaluating pathogen whole genome sequencing in public health

High-precision and cost-efficient sequencing for real-time COVID-19 surveillance

Genomic characterisation and epidemiology of 2019 novel coronavirus: Implications for virus origins and receptor binding

A pneumonia outbreak associated with a new coronavirus of probable bat origin

A comparison of whole genome sequencing of SARS-CoV-2 using amplicon-based sequencing, random hexamers, and bait capture

COVseq is a cost-effective workflow for mass-scale SARS-CoV-2 genomic surveillance

Multiple approaches for massively parallel sequencing of SARS-CoV-2 genomes directly from clinical samples

Next generation sequencing of SARS-CoV-2 genomes: challenges, applications and opportunities

Shotgun metagenomics, from sampling to analysis

Oligonucleotide capture sequencing of the SARS-CoV-2 genome and subgenomic fragments from COVID-19 individuals

A rapid, cost-effective tailed amplicon method for sequencing SARS-CoV-2. BMC Genom

Touchdown" PCR to circumvent spurious priming during gene amplification

High and low annealing temperatures increase both specificity and yield in touchdown and stepdown PCR

Selected insights from application of whole-genome sequencing for outbreak investigations

The next phase of SARS-CoV-2 surveillance: Real-time molecular epidemiology

Genomic analysis of viral outbreaks

The impact of real-time whole genome sequencing in controlling healthcare-associated SARS-CoV-2 outbreaks

Characterization of SARS-CoV-2 infection clusters based on integrated genomic surveillance, outbreak analysis and contact tracing in an urban setting

Evaluation of NGS-based approaches for SARS-CoV-2 whole genome characterisation

Rapid genomic characterization of SARS-CoV-2 viruses from clinical specimens using nanopore sequencing

Rapid high-throughput whole-genome sequencing of SARS-CoV-2 by using one-step reverse transcription-PCR amplification with an integrated microfluidic system and next-generation sequencing

Whole genome sequencing of SARS-CoV-2: Assessment of the ion torrent ampliseq panel and comparison with the illumina-MiSeq ARTIC protocol

Real-time SARS-CoV-2 diagnostic and variants tracking over multiple candidates using nanopore DNA sequencing

High throughput detection and genetic epidemiology of SARS-CoV-2 using COVIDSeq next-generation sequencing

Mutations in emerging variant of concern lineages disrupt genomic sequencing of SARS-CoV-2 clinical specimens

Improvements to the ARTIC multiplex PCR method for SARS-CoV-2 genome sequencing using nanopore

Disentangling primer interactions improves SARS-CoV-2 genome sequencing by multiplex tiling PCR

Fast and accurate short read alignment with Burrows-Wheeler transform

An amplicon-based sequencing framework for accurately measuring intrahost virus diversity using Pri-malSeq and iVar

The Sequence Alignment/Map formt and SAMtools

Designed the study, participated in viral genome sequencing, assay development and validation, data analysis and wrote the manuscript. A.Ge. designed, developed and validated the bioinformatics analysis pipeline

participated in study design and validation, viral genome sequencing, sample collection, assay validation and preparation and editing of the manuscript. K.L. Designed, developed and validated the high-throughput laboratory automation used for viral genome sequencing in this study. R.M.K designed and developed the viral genome database used in this study, implemented automated lineage and clade assignments to populate the database, performed data analysis and assisted with sample collection and edited the manuscript. B.A. assisted with high-throughput laboratory automation design and development and participated in viral genome sequencing

Designed, developed and deployed the informatics infrastructure and software necessary to operationalize the assay pipeline and to provide and interface to the laboratory scientists who performed the testing and reporting. P.T. and F.L. provided medical oversight and assisted in securing necessary resources and funding for this study

NovaSeq runs and sequence IDs are listed in Supplementary Materials file "Sequenced_sample_list.xlsx". Fasta consensus sequences in compressed (zip) file format will be provided upon request. Large scale surveillance sequencing was performed for the CDC National SARS-CoV-2 Strain Surveillance (NS3) project (https:// www. aphl. org/ progr ams/ prepa redne ss/ Crisis-Manag ement/ COVID-19-Respo nse/ Pages/ Seque nce-Based-Surve illan ce-Submi ssion. aspx). Consensus sequences were submitted to the CDC and were subsequently deposited by the CDC in the GISAID database (https:// www. gisaid. org/); originating lab: Quest Diagnostics Incorporated. Sequences submitted by Quest Diagnostics as part of this project may be retrieved from the GISAID search interface by using a virus name search, keyword "QDX".

This study was funded by Quest Diagnostics.

All authors were employees of Quest Diagnostics during the time this study was being conducted. Quest Diagnostics is a reference laboratory that provides commercial testing for SARS-CoV-2, including whole genome sequencing.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-022-06091-0.Correspondence and requests for materials should be addressed to R.M.K. or F.L.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.