{'fraction_non_alphanumeric': 0.04252278006074683, 'fraction_numerical': 0.03229408611756298, 'mean_word_length': 4.708567057623662, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of SARS-CoV-2 Variants of Concern Using Amplicon Next-Generation Sequencing

ABSTRACT COVID-19 is caused by SARS-CoV-2, several virulent variants of which have emerged since 2019. More than 529 million people have been infected, and at least 6 million have died. Our aim was to develop a fast, accurate, low-cost method for detecting and identifying newly emerging variants of concern (VOCs) that could pose a global threat. The 341-bp DNA sequence of a specific region of the SARS-CoV-2’s spike protein was amplified by a one-step PCR on RNA samples from 46 patients. The product was sequenced using next-generation sequencing (NGS). DNA sequences from seven genomes, the original Wuhan isolate and six different representative variants obtained from the GISAID website, were used as references. Complete whole-genome sequences from local isolates were also obtained from the GISAID website, and their RNA was used for comparison. We used an amplicon-based NGS method (termed VOC-NGS) for genotyping and successfully identified all 46 samples. Fifteen (32.6%) were like the original isolate. Twenty-seven were VOCs: nine (19.5%) Alpha, eight (19%) Delta, six (14%) Beta, and four (8.7%) Omicron. Two were variants of interest (VOI): one (2%) Kappa and one (2%) Zeta. Two samples were mixtures of two variants, one of Alpha and Beta and one of Alpha and Delta. The Spearman correlation between whole-genome sequencing (WGS) and VOC-NGS was significant (P < 0.001) with perfect agreement (Kappa = 0.916) for 36/38 (94.7%) samples with VOC-NGS detecting all the known VOCs. Genotyping by VOC-NGS enables rapid screening of high-throughput clinical samples that includes the identification of VOCs and mixtures of variants, at lower cost than WGS. IMPORTANCE The manuscript described SARS-Cov-2 genotyping by VOC-NGS, which presents an ideal balance of accuracy, rapidity, and cost for detecting and globally tracking VOCs and some VOI of SARS-CoV-2. A large number of clinical samples can be tested together. Rapid introduction of new mutations at a specific site of the spike protein necessitates efficient strain detection and identification to enable choice of treatment and the application of vaccination, as well as planning public health policy.

B.1.427 and B.1.429 lineages, and in B.1.617 (Kappa) lineages and sublineages, and E484K, seen in all B.1.525 (Eta) and P.2 (Zeta) variants but only in some strains of the B.1.526 (Iota) variant. Based on the CDC data available as of preparing the manuscript, there are no variants of SARS-CoV-2 that fall in the VOHC category (6,17). In this study, VBMs and VOCs will be designated "VOC." A quick and accurate but low-cost surveillance system for detecting, identifying, tracking, and controlling the emergence of VOC variants is needed. Real-time identification of VOC and VOI will have a significant impact on the management of the current pandemic and on patient care (e.g., monoclonal antibody therapy) (24). Currently, applying real-time identification requires WGS capacity, which might not exist or might be restricted with regard to high-throughput capacity. WGS is also expensive and time-consuming. An accurate amplicon-specific target screening strategy would allow full screening for each positive sample. Several assays have been employed to detect only specific VOCs. For example, Aoki et al. (25) developed a genotyping platform for SARS-CoV-2 variants using high-resolution melting (HRM) analysis. Reverse transcription-quantitative PCR (RT-qPCR) assays, incorporating 14 primers and fluorescence-labeled probes, have also been developed for rapid detection of the SARS-CoV-2 mutations N501Y, 69 to 70del, K417N, and E484K, which occur in variants having a clinical impact (26). However, these technologies are specifically designed for previously identified variants. Here, we developed a quick, low-cost amplicon-based NGS method (termed VOC-NGS) that can maximize the testing of variants and broaden clinical laboratories' ability to actively participate in identifying new VOCs.

RESULTS
The PCR amplification produced two fragments with the combination of three primers (Table 1), two forward and one reverse, and was selected to amplify the two nucleotide sequences seen as the bands 213 bp and 421 bp (Fig. 1). A lower concentra- tion (0.125 mM compared to 0.5 mM) of the second forward primer, HS478F, was used to control excessive amplification of the lower-molecular-weight (213 bp) product, at the expense of the higher-molecular-weight product (421 bp). VOC-NGS experiments were performed on RNA samples from 46 local patients, 10 of which were selected as they were known to contain the Alpha variant and the original strain. In parallel, and for comparison, WGS was performed. The WGSs of another 10 samples ( Table 2) that were already in the GISAID website were used for comparison in assessing the efficacy of genotyping by VOC-NGS introduced here, which was able to identify all 46 (100%) samples studied. The number of variants found among the 46 samples studied were as follows: 15 (32.6%) were like the original isolate, 9 (19.5%) Alpha, 8 (19%) Delta, 6 (14%) Beta, and 4 (8.7%) Omicron. In addition, two VOIs were revealed, one (2.4%) Epsilon and one (2.4%) Zeta. Interestingly, two samples showed a mixture of different variants, one a mixture of Alpha and Beta variants (sample 824), the other a mixture of Alpha and Delta variants (sample 845). These mixtures were confirmed by WGS. The association of WGS and VOC-NGS was shown to be highly significant by Spearman correlation P , 0.001, which was applied to the 38 samples evaluated by both assays where 36 of 38 (94.7%) perfect agreements (Kappa = 0.916) occurred. There was 100% identity between the results of WGS and the genotyping by our VOC-NGS amplicon-based method. This method can detect any new mutation within this specific HSR. It can also reveal other VOIs such as Gamma (K417T, E484K, N501Y), Epsilon (L452R), and Zeta (E484K) ( Fig. 2 and Table 2).
Furthermore, genotyping by VOC-NGS was able to detect and quantify coinfection-contamination of different variants as shown by the clinical samples 824 (Alpha, 35%; Beta, 65%) and 845 (Alpha, 47%; Delta, 53%) in Fig. 3 where the percentages indicate the degree of coinfection-contamination. In Fig. 3, at position 283, sample 824 has the mutation N501Y with the nucleotide A changes to T at 98% coinfection contamination, and at position 232, sample 845 has the mutation E484K with the nucleotide G changed to A at 65% coinfection contamination (see the arrows below the chart and the box for quantification of the conversion at the top of the chart). The chart was displayed with local integrative genomics viewer (IGV).

DISCUSSION
NGS was used in developing a new assay, VOC-NGS genotyping, for use in standard clinical laboratories to rapidly, reliably, and inexpensively screen for VOCs of SARS-CoV-2. This method is able to detect and identify the different spike protein mutations that have occurred in VOCs with a 94.1% agreement with identifications of local variants done by WGS. With a kappa value of 0.916, this is considered to be almost perfect agreement since a kappa value of between 0.81 and 1.00 is accepted as almost perfect agreement (27). VOC-NGS genotyping successfully identified the original isolates and the Alpha, Beta, Delta, and Omicron variants from the four consecutive surges of the COVID-19 pandemic in the region.
Several assays have been employed to detect and identify only VOCs, e.g., high-resolution melting (HRM) analysis to detect and identify variants of SARS-CoV-2 (25) and RT-qPCR assays to detect the mutations N501Y, 69 to 70del, K417N, and E484K occurring in variants of SARS-CoV-2 that have a clinical impact, using 14 sets of primers and fluorescence-labeled probes (26). Differences in annealing temperature can affect both assays. Both assays require a positive control for detecting the presence of each mutation and for building different HRM graphs. Also, both assays provide only limited information with regard to specific probes, e.g., the targeting of the N501Y mutation by its specific probes. Changes close to this mutation, a region known to be polymorphic, could prevent the probe's annealing to the target site and consequently miss detection of this variant. The same goes for the HRM assay, where any change in the annealing temperature of the product component would cause it to change its melt curve and produce misleading peaks. Furthermore, appropriate positive controls will be needed for each SNP in each run. In contrast, genotyping by VOC-NGS is unbiased and  provides entire sequences despite the occurrence of mutations within the product. If a mutation occurs at the target site where the primers anneal, a product will not be produced. Such events could be subject to follow-up investigation by longer PCR and Sanger sequencing. As genotyping by VOC-NGS is not based on specific probes and the need of doing WGS, it should enable early global tracking and monitoring of the spread of the known and new introduced coronavirus variants. Owing to its rapidity, accuracy, and high throughput, its wide employment would alert people in authority, policy makers and public health providers, to the emergence, spread, and reemergence of the virus and its variants, enabling them to try to prevent their spread.   There are advantages to using VOC-NGS genotyping, an amplicon-based NGS method, over WGS using the ARTIC v3 system (28). There is no need for amplification of the whole genome using two PCR pools, followed by Nextera DNA Flex preparation, which increases costs (ARTIC v3 sequencing costs £975 per 100 samples, using CoronaHiT-Illumina sequencing of the about 30,000 bp of the SARS-CoV-2 genome) (29). The total cost for genotyping by VOC-NGS is only £180 for 100 samples, which includes sequencing at a depth of 400Â with a 150 MID kit costing £30, Luna one-step ready-mix costing £130, and Agencourt AMPure XP and Nextera barcodes costing £20. In total, this amounts to almost a fifth (18%) of the cost of Corona WGS (£180 versus £975).
While WGS provides more information on the virus and its variants, these fiscal calculations clearly show that detection and identification of variants of COVID-19 using genotyping by VOC-NGS costs less than that by WGS. Also, genotyping by VOC-NGS, which can be operated in different setups, is fast, as it produces cDNA and carries out specific HSR amplification with the specific primers in a single reaction. Moreover, free online analysis tools such as Galaxy, Genome Detective, and IGV produce accurate consensus sequences of the 341-bp region and can be employed by any standard laboratory performing data analysis, thus avoiding bioinformatic analysis costs. The products prepared from the sample for genotyping by VOC-NGS are suitable for use in all Illumina platform machines (MiniSeq, MiSeq, Nextseq, Iseq100i, and Novaseq machines) for 150-bp single and paired-end reads. For high-throughput screening, 384 samples can be loaded and run together owing to the barcode combinations of indices that can be multiplexed together, and the protocol is compatible with a robotic system.
Genotyping by VOC-NGS was able to identify all 341 bp of the spike HSR. This is where all the mutations giving rise to the known VOCs occur and encompasses the var- However, the variants VOI Theta GR/1092K.V1 (P.3), first detected in the Philippines, and VOI Iota GH/253G.V1 (B.1.526), first detected in New York, USA, cannot be detected with genotyping by VOC-NGS, as their underlying mutations fall outside the HSR target zone. In this case, new PCR primers can be developed to detect and identify the mutations of new VOIs situated outside the HSR that can be included in the same PCR. Thereafter, the genotyping procedure would remain as described for VOCs. The cost of preparing these additional primers would be very low. The combination of primers used in the method described in this article is more straightforward in detecting and identifying VOCs and most VOIs and, therefore, has greater diagnostic utility than single-mutation-specific PCR approaches like the RT-qPCR and HRM assays.
The VOCs affect the distribution, expansion, and reemergence of the COVID-19 pandemic during its different surges and its disease severity and vaccine efficacy (12,20,21,30). Genotyping by VOC-NGS is also able to detect and identify the Delta plus variant of SARS-CoV-2 recently discovered in Israel since the mutations K417T, E484K, and L452R are in the target region of the variant's spike protein. VOC-NGS genotyping was able to reveal and quantify a mixture of two-variant infection as shown in two samples (824 and 845) in this study. This was previously reported in a rare case of coinfection of an unvaccinated elderly woman who was shown to be infected at the same time with the Alpha and Beta variants of SARS-CoV-2 (31). Furthermore, recent study on SARS-CoV-2 strains reveals potential contribution of coinfection with and recombination between different strains to the emergence of new strains (32).
To check for new mutations, the sampling of screen-negative specimens for WGS will be an important component of a full surveillance strategy, as the approach described here cannot predict the existence of new variants that carry new mutations. Addition of more amplicons by targeting new regions of interest for use in genotyping by VOC-NGS is a good option with no significant price difference.
This method could play a role in choosing therapy as an alternative to administering bamlanivimab and etesevimab, which are generally used to treat cases of COVID-19 but have shown reduced efficacy against the Beta variant B.1.351 that contains the mutations K417N plus E484K plus N501Y (23). Identifying these mutations will indicate the need for alternative therapy for patients infected with this variant (23,30).
In cases of low RNA abundance (high RT-PCR cycles), the sensitivity of VOC-NGS is higher than that of NGS. This was apparent in clinical sample S747. This sample had threshold cycle levels (32), and failed in WGS, but was clearly identified by VOC-NGS genotyping. This is due to the targeting specific amplicons on the genome.
The technique's main limitation remains the need of an NGS facility. However, this can be solved by shipping DNA libraries to either local or international companies, such as Macrogen, https://www.macrogen-europe.com/, and Novogene, https://en.novogene.com/, or laboratories with NGS capabilities.
In summary, genotyping by VOC-NGS presents an ideal balance of accuracy, rapidity, and cost for detecting and globally tracking VOCs and some VOIs of SARS-CoV-2. A large number of clinical samples can be tested together. Rapid introduction of new mutations at a specific site of the spike protein necessitates efficient strain detection and identification to enable choice of treatment and the application of vaccination, as well as planning public health policy.

MATERIALS AND METHODS
RNA samples. This study was part of the officially permitted diagnostic procedure of The Hebrew University Hadassah Medical Center; therefore, no additional approval from the Institutional Review Board was required. Samples were collected in Jerusalem between February and July 2021. Table 2 describes the clinical RNA samples and viral variants used in this study. The genetic sequence of the spike protein of SARS-CoV-2 was that of the Wuhan-Hu-1 isolate (GenBank accession number: NC_045512.2, 21563 to 25384). The nucleic acid sequence between the mutations K417N and N501Y was targeted as it encompasses all the mutations that have produced VOCs (Table 2 and Fig. 4). The schematic overview of workflow is shown in Fig. 5.
Whole-genome sequencing of viruses from clinical samples. RNA was extracted using the MagDEA Dx SV kit on the MagLEAD platform (Precision System Science Co., Ltd.) according to the manufacturer's instructions. A volume of 280 mL lysis buffer was added to a 220-mL sample and eluted in 50 mL. The RNA concentration/quality was determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher), and extracted RNA was kept at 270°C for further analysis.
Sequencing libraries were prepared using the Illumina COVIDSeq test and sequenced on Novaseq 6000 producing at least 3.3 million paired-end reads (150 nucleotides [nt]) per library. Library quality was analyzed using FastQC (version 0.11.8, Babraham Bioinformatics). Reads were aligned to the genome using Bowtie2 (version 2.3.4.3) with the command options -k 4 -no-discordant. Reads with more than 6 variants in 100 bases were discarded (single nucleotide variants [SNV], deletion or insertion each count as one variant). Variants were called using ivar variants (version 1.3.1). Consensus sequence was built based on the ivar variants table using the R Biostrings package according to the following rules. Positions with fewer than 10 reads were called N. Variants with frequency higher than 0.7 were included in the consensus sequence. Variants with frequency between 0.3 and 0.7 and at least 50 reads were considered "wobbles" using the IUPAC letters. Consensus sequences were submitted to Pangolin command-line tool (pangolin version 2.2.2 and pangoLEARN version 2021-02-12) and Nextclade (version 0.12.0) to determine the PANGO lineage and clade.
PCR amplification of the HSR region. Two sets of primers were designed, using the primer 3 online program, and synthesized. In accordance with the Nextera Illumina DNA library preparation kit (Illumina, San Diego, CA, USA), they included specific Illumina adaptor overhang sequences added at the 59 end: TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG for the forward primer and GTCTCGTGGGCTCGGAGATGTGTA TAAGAGACAG for the reverse primer (Fig. 5).
PCR amplification to obtain the two required fragments was done with a combination of three primers (Table 1), two forward primers, MutF1 and HSF487, and one reverse primer, MutR1. This resulted in amplification of two nucleotide sequences (Fig. 1) detected as 213 bp (a 146-bp sequence flanked by 67 bp of primers HSF487 and MutR1) and 421 bp (a 354-bp sequence flanked by 67 bp of primers MutF1 and MutR1). To avoid an excess of the lower-molecular-weight product, a lower concentration of the second forward primer was used (0.125 mM HSF487 as opposed to 0.5 mM MutF1 and MutR1).
We conducted reverse transcription and PCR amplification as a one-step reaction that comprised cDNA preparation and PCR amplification, using the Luna Universal Probe one-step RT-qPCR kit (E3006, NEB) according to the manufacturer's protocol. A total of 2.5 mL extracted RNA in a reaction volume of 20 mL was used. The PCR conditions were a reverse transcription cycle of 10 min at 55°C, then 1 min at 95°C, and then 40 cycles of 10 s at 95°C, together with annealing and extension in one step of 30 s at 60°C. PCR products were purified, using AMPureXP beads (Beckman Coulter, Brea, CA, USA). Twenty-five microliters of beads was added to the PCR product, incubated for 5 min at room temperature, and then put on a magnetic plate for 2 min, followed by two washes with freshly made 80% ethanol and then drying for 2 min on the plate. Finally, beads were eluted with 22 mL of elution buffer. The purified amplicons in a volume of 10 mL underwent eight cycles of PCRs to anneal dual-indexed barcodes with unique sequences: index 1 read 59-CAAGCAGAAGACGGCATACGAGAT[i7]GTCTCGTGGGCTCGG, index 2 read 59-AATGATACGGCGACCACCGAGATCTACAC[i5]TCGTCGGCAGCGTC. i7 and i5 are unique 8-bp sequences for identifying amplicons derived from each sample (Illumina, 2013). Then, 10 mL of each sample was pooled and purified again with AMPureXP beads. The purified pool product was adjusted to 4 nM and sequenced with a NextSeq 500/550 mid output kit v2.5 (150 cycles), using a NextSeq 500 machine (Illumina, San Diego, CA, USA), aiming at ;1,000 reads giving a depth of 400Â per sample. Binary base call (BCL) output files from a Nextseq 500 machine were converted to FASTQ format, using BCL to FASTQ (bcl2fastq v2.20.0.422, Illumina, Inc.). FASTQ format files were analyzed using the Galaxy program (Galaxy version 0.7.17.1) (https://usegalaxy.eu/) (33). The 341-bp DNA sequence was obtained from the sequencing of two PCR products. The first 150-bp DNA sequence was from read 1 (forward primer 1) from PCR1, followed by a gap of 40 Ns, and then the second 150-bp sequence started from forward Corona VOC Detection and Identifications by NGS Microbiology Spectrum primer 2 (read 1) from PCR2 (see Fig. 2). The 341-bp consensus sequences were obtained by mapping reads with BWA-MEM map medium and long reads (.100 bp) against the Wuhan-Hu-1 SARS-CoV-2 reference isolate (selected 341 bp) fragment from isolate hCoV-19/Wuhan/Hu-1/2019, GenBank accession number NC_045512, GISAID accession ID EPI_ISL_402125. The mapping was displayed with a local integrative genomics viewer (IGV), and the consensus sequence from the IGV was copied and accepted as the sequence of the sample's variant. The genome detective virus tool version 1.133 (https://www .genomedetective.com/app/typingtool/virus/) was used to directly reveal the spike mutation at the DNA and amino acid levels (34). A free online multiple sequence alignment program (35) was used to align these HSR DNA sequences of the spike protein from patients with COVID-19 and different VOCs from the GISAID website ( Fig. 4 and Table 2). Statistical analysis. The Statistical Package for the Social Sciences (SPSS) computer software program version 15.0 (SPSS Inc.) was used to test the Spearman correlation as a nonparametric test used to measure the degree of association between VOC-NGS assay and the WGS assay results. The free online GraphPad (https://www.graphpad.com/quickcalcs/kappa1/?K=5) program was used for kappa agreement tests (27). Cohen's kappa coefficient (k ) is a measure of the agreement between two tests beyond that expected by chance, where 0 is chance agreement and 1 is perfect agreement for the VOC-NGS assay performance against the WGS method.
Data availability. The original data using a next-generation sequencing data set from the GISAID website used in the study are found in Table 2. The data sets used and/or analyzed during the present study are available from the corresponding author on reasonable request.