Validation of a new strategy for the identification of SARS-CoV-2 variants by sequencing the spike gene by Sanger

Introduction The emergence of multiple variants of SARS-CoV-2 during the COVID-19 pandemic is of great world concern. Until now, their analysis has mainly focused on next-generation sequencing. However, this technique is expensive and requires sophisticated equipment, long processing times, and highly qualified technical personnel with experience in bioinformatics. To contribute to the analysis of variants of interest and variants of concern, increase the diagnostic capacity, and process samples to carry out genomic surveillance, we propose a quick and easy methodology to apply, based on Sanger sequencing of 3 gene fragments that code for protein spike. Methods Fifteen positive samples for SARS-CoV-2 with a cycle threshold below 25 were sequenced by Sanger and next-generation sequencing methodologies. The data obtained were analyzed on the Nextstrain and PANGO Lineages platforms. Results Both methodologies allowed the identification of the variants of interest reported by the WHO. Two samples were identified as Alpha, 3 Gamma, one Delta, 3 Mu, one Omicron, and 5 strains were close to the initial Wuhan-Hu-1 virus isolate. According to in silico analysis, key mutations can also be detected to identify and classify other variants not evaluated in the study. Conclusion The different SARS-CoV-2 lineages of interest and concern are classified quickly, agilely, and reliably with the Sanger sequencing methodology.


Introduction
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) Q4 is a coronavirus that emerged in 2019 and is responsible for the current pandemic known as coronavirus disease 2019 (COVID-19) 1 . COVID-19 has since become a global public health threat, causing approximately 405 million infections and 5.8 million deaths worldwide 2 .
Early diagnosis is crucial to controlling the spread of COVID-19. The molecular detection of SARS-CoV-2 nucleic acid by real-time polymerase chain reaction (PCR) is considered the gold standard. However, it can be affected by virus mutations in the 5 (UTR), ORF1ab, spike (S), membrane, and nucleocapsid regions 3,4 , which favour an increase in the rate of false negatives 5 .
Mutations in the SARS-CoV-2 genome generate changes in protein structures, affecting their virulence, transmissibility and antigenicity 6,7 . Therefore, genomic surveillance of the virus has been essential in detecting early variants of interest (VOI) or variants of concern (VOC). All variants have mutations in the protein S 4 gene. Therefore, most therapeutic strategies and vaccine development have focused on this gene.
The detection of the different mutations and the typification of the variants is carried out by next-generation sequencing (NGS) of the complete SARS-CoV-2 genome. The virus sequences are deposited in the Global Initiative on Sharing All Influenza Data EpiCoV TM database, so that genetic information from positive samples for SARS-CoV-2 8 can be shared. Similarly, there are platforms that enable genomic and phylogenetic analyses of the variants, such as Nextstrain 9 .
Although NGS techniques are crucial for molecular surveillance, they require a high degree of technical expertise, knowledge of bioinformatics and access to expensive equipment, rendering widespread use difficult. For these reasons, other methodologies have been designed, such as Sanger sequencing and real-time PCR, as variant screening tests 3,10 .
In Colombia, SARS-CoV-2 genomic characterisation and surveillance activities are centralised in the National Institute of Health. However, with the aim of contributing to genomic surveillance, increasing diagnostic capacity and simple workflow, and improving variant tracking in the shortest possible time, we propose a method based on Sanger sequencing of three fragments that partially cover the gene encoding protein S to detect VOI and VOC.

In silico simulation of detection of variants of interest and concern
The possible mutations in the amino acid sequence of the Gen-Bank protein S were simulated. The reference sequences of the VOC (alpha, beta, gamma, delta and omicron) and VOI (lambda and mu) were obtained. They were aligned with the Sequencher 5.4.6 program (Gene Codes Corporation) and the variants were identified with the Nextstrain-Nextclade Web 1.10.0 web-based program, Nextclade CLI 1.7.0 (09-12-2021) (https://clades.nextstrain.org/). RNA extraction and cDNA synthesis RNA was obtained using automated extraction systems with magnet technology, KingFisher TM Duo Prime System (Thermo Fisher Scientific Inc., USA.) with the MagMAX TM Viral/Pathogen Nucleic Isolation kit (Thermo Fisher Scientific Inc., USA). Similarly, manual extractions of some samples were carried out with the Quick-RNA Viral Kit (Zymo Research, USA). Total RNA was converted to cDNA using the TruScript TM First Strand cDNA Synthesis kit (Norgen, Thorold, Canada), following the manufacturer's recommendations.

Sanger sequencing and variant analysis
The PCR products were quantified, diluted to 40 ng/l and purified with ExoSAP-IT TM for PCR Product Clean-Up (Applied Biosystems, Lithuania). Subsequently, they were sequenced in both directions using the BigDye ® Terminator v3.1 Cycle Sequencing kit (Applied Biosystems/Hitachi High Technologies, Japan) and purified with the BigDye ® XTerminator TM Purification Kit. Finally, capillary electrophoresis was performed with the 3500 Genetic  Analyzer equipment (Applied Biosystems, Lithuania). The consensus sequences were generated using Sequencher 5.4.6 (Gene Codes Corporation, USA).
To confirm the results obtained from the Sanger sequencing, the samples were also analysed at the Sequencing Unit of the University of Antioquia, with support from the PECET (Programa de Estudio y Control de Enfermedades Tropicales [Programme for the Study and Control of Tropical Diseases)]-Immunovirology groups of the same university, by means of NGS using the iSeq TM 100 (Illumina, USA), MiSeq TM (Illumina, USA) and Oxford Nanopore MinION (Oxford Nanopore Technologies, United Kingdom) systems. Variant nomenclature was performed according to the Nextstrain-Nextclade Web 1.10.0 online program, Nextclade CLI 1.7.0 (2021-12-09) and PANGO Lineages (https://cov-lineages.org/).

Statistical analyses
The mutation data were tabulated and analysed using the GraphPad Prism v7.04 program. The concordance of the determination of variants between Sanger sequencing and by NGS was calculated using the kappa index, considering NGS as the gold standard and defining two possible categories (present or absent) in relation to the mutations or changes evidenced in the S protein.

Results
Three pairs of primers were designed that start at position 21,620 and end at position 25,076 of the reference sequence NC 045512 of the S gene of the SARS-CoV-2 virus, which are listed in Table 1. The amplicons generated overlapped each other and generated a product of 3434 bp with a theoretical coverage of 89.9% of the S gene (3822 bp), as can also be seen in Fig. 1.
According to the in silico analysis, the mutations present in the variants reported to date do not interfere with the yield and amplification of the fragments. However, a single nucleotide change was identified in the LIMEF3 FW primer that amplifies fragment 3.

In silico analysis of variants of interest and concern
In silico analysis of the VOI and VOC reported to date in the World Health Organization's Tracking SARS-CoV-2 variants was conducted. Theoretically, with the Sanger methodology, the specific mutations S:G75V, S:T76I, S:F490S and S:T859N, and seven deletions of the 21G variant (lambda) corresponding to the C.37 lineage can be detected. In the case of the mu variant, the relevant S:R346K and S:Y144S mutations can be detected, as well as the seven mutations also reported in the other VOI and VOC.
For the VOC, in alpha, the Sanger results identified the specific mutations S:A570D, S:T716I and S:S982A. However, the mutation located at amino acid 1118 (S:D1118H) could not be detected. In the case of the beta variant, mutations S:D80A and S:D215G and deletions S:L241-, S:L242-and S:A243-can be detected, and in gamma, variants S:D138Y, S:R190S and S:T1027I can be detected. For the delta variant, the mutations shared by clade 21A and 21J (S:R158G, E156-and S:F157-) could be detected, but subclade 21A could not be differentiated from 21J due to the presence of mutations in the genes of the nucleocapsid and ORF1a-b, which are not analysed by the Sanger methodology. Finally, in the case of the omicron variant, two mutational profiles close to each other were found, clades 21K and 21L. For the 21K clade, eight specific mutations could be detected, but in the 21L clade only seven mutations would be detected, which would allow differentiation between the two clades.

Analysis of variants in positive samples for SARS-CoV-2
The amplicons obtained from the three conventional PCRs experimentally cover 3434 nucleotides of the S gene, corresponding to 1144 amino acids. Good quality Sanger sequences were obtained in all 15 clinical samples. Two alpha samples, three gamma, one delta, three mu and one omicron sample were classified and five strains were related to the initial isolate of the Wuhan-Hu-1 virus.
In the blinded comparison of the results of the Sanger sequencing and NGS (gold standard), the classification results of both methods in the 15 samples were concordant (Appendix B Supplementary material, annex 1).
Samples B1 (OM189518) and 5703 (OM128430) were classified in clade 19B, while the NGS result for both samples placed them in 20A. The mutations found were S:D614G for sample B1, and the mutations W152R and D614G were found in 5703. For their part, the sequences with codes 6041 (OM060697), 2403 (OM189474) and 12,422 (OM127616) were diagnosed in early 2021 and classified with pangolineage B.1.625, clade 20A.
Sample 36,211 (OM142697) was classified by Sanger as a delta variant, but was located in clade 21A, while it was located in 21J by the NGS methodology. Both sequencing methods present the same S protein mutation profile. Sample M1 (OM066778) was categorised as an omicron variant in the 21K clade, with more than 36 mutations detected in the S protein amino acid sequence.

Discussion
The comparison of the results of the Sanger and NGS sequencing of the COVID-19 samples evaluated showed that the Sanger methodology achieves the correct classification of different lin- eages of SARS-CoV-2 with results consistent with those obtained by NGS (Supplementary material, annex 1).
The statistical analysis with the kappa index was 0.146 (95% CI −0.032 to 0.325), showing a slight agreement between methods, probably associated with the fact that none of the methodologies detect all the mutations reported in the Tracking SARS-CoV-2 variants database. It should be noted that the Sanger methodology detects key mutations to identify the clade to which a strain belongs.
Although technically it is not possible to amplify and sequence the 3822 bases of the S gene using the Sanger methodology in a sin-gle reaction 13 , in our study we proposed the amplification of three overlapping fragments to achieve 90% coverage of the S gene. This is efficient compared to other sequencing methodologies 13,14 , since it shortens processing time, does not require a minimum number of samples for analysis, reduces the number of reactions necessary to obtain results and reduces the associated costs.
The methodology proposed in this article requires adequate collection, conservation and transport of the samples, which ensures a complete RNA with cycle threshold values less than 25. These data are similar to those recommended in other investigations 13,15 . The protocol was tested with samples extracted by automated methods of magnetic beads and extractions with silica columns, generating good results, particularly with the column methodology.
Although not all the SARS-CoV-2 variants were found in the clinical samples, the in silico analysis of protein S showed that this methodology guarantees the detection of the amino acids of interest S:K417, S:L452, S:S477, S:T478, S:E484, S:F490, S:S494 and S:N501. These mutations are of great importance because the WHO and the United States' CDC defined them as key mutations present in VOI and VOC for genomic surveillance purposes 13 .
The S:D614G mutation was found in all the samples analysed. It is called the universal mutation and all VOI and VOC are considered to be derived from it 16 . S:D614G is related to increased viral replication in the upper respiratory tract and increases the infection's efficiency 15 . This adaptive advantage is related to the increase in COVID-19 cases in some European countries in mid-2020 and the subsequent worldwide spread 17,18 .
The alpha variant is associated with multiple mutations in the S protein, the most prominent one being S:N501Y and the 69-70 deletion. The latter generates a failure of some real-time PCR assays 5,19 and is found in other variants, such as 20A, 21K (omicron) and 21L (omicron).
Based on covariate data, mutations S:A570D, S:T716I, S:S982A and S:D1118H are associated with the alpha variant 16 . In this study, the S:D1118H mutation could not be detected due to lack of coverage in that region. However, the mutations identified are sufficient to categorise the variant.
Gamma is associated with multiple mutations; the specific ones are S:N501Y, S:E484K, S:L18F, S:K417T and S:H655Y. They are shared by several VOC, except the delta variant 16 . Our methodology is able to detect around eight of the characteristic mutations, including three specific mutations: S:D138Y, S:R190S and S:T1027I.
The delta variant and the omicron variant are considered to be VOC, and their rapid spread has been the cause of the outbreaks since mid-2021 20 . The delta variant has two subclades; the Sanger and NGS results analysed in Nextstrain generated the classification 21A and 21J, respectively, for sample 36211. The discrepancy obtained by the Sanger and NGS methods is due to the fact that the 21J subclade contains all the mutations of the 20A clade with additional mutations in the N gene (N:G215C), as well as in ORF1a (ORF1a:A1306S, ORF1a:V2930L, ORF1a:T3255I, ORF1a:T3646A) and ORF1b (ORF1b:A1918V, and ORF7b:T40 I) 16,21 , which are not analysed in this study.
The omicron variant emerged in November 2021 and early infections and sequence data are predominantly from South Africa 22 . The identification of this variant is very significant because it is currently the most common one, its mutations are associated with greater transmissibility, immune evasion and greater binding affinity to ACE2 23,24 . It shares two mutations with the alpha VOC but has around 26 specific mutations 16 , which are mostly detected by the Sanger method.
The mu VOI emerged in early 2021 and was very common in Colombia 16 . This variant has mutations that are also found in VOC and VOI such as S:T95I, S:E484K, S:N501Y, S:P681H, and S:D95N. The S:R346K mutation has been considered specific for the mu variant, although this mutation also occurs in the omicron variant.
It was not possible to obtain samples from the beta and lambda variants to perform the experimental analysis; the in silico analyses confirm that the Sanger sequencing method can detect characteristic mutations.
Our group acknowledges that the number of samples analysed is low, but it allows us to demonstrate that the methodology has the capacity to experimentally categorise VOC and VOI using freely accessible online applications without requiring the creation of subjective algorithms. For further studies, we recommend a larger number of samples and other variants not included in the study.
This methodology overcomes some of the limitations of 25 , intercalating 26 and allele-specific real-time probe PCR methods 27 , where the use of specific primers and/or probes is necessary to evaluate each mutation and with prior knowledge of circulating mutations. Additionally, some of these PCR are not specific for determining genetically related variants (subclades), they do not allow the appearance of new variants to be evaluated and they underestimate the real number of circulating variants.
Recently, Sanger sequencing methodologies that study the S protein have been published and even standardised commercial protocols are already available [28][29][30] . Several of those studies require a larger number of primers targeting domains of interest, such as receptor-binding domains and receptor-binding motifs 30 , as well as analysing mutations in other regions 28 , and in some cases they perform a total sequencing of the S protein. All these previous methodologies involve performing multiple PCRs 29 , which increases the cost compared to the methodology that we propose here.
It is important to stress that this approach does not replace NGS technology and other PCR-based assays. Rather, it seeks to increase the capacity of the epidemiological surveillance network, being used as a variant screening method, and in cases in which NGS has little coverage or depth it can be used as a strategy to complete the sequences.
In conclusion, this methodology enables the detection of mutations in the S protein in a streamlined and precise way, quickly and with infrastructure that is available in clinical, public health and research laboratories in Colombia that serve as support for the genomic diagnostic and surveillance network of the SARS-CoV-2 virus.

Funding
This work was financed with resources from Colombia's general awards system (BPIN 2020000100152) and from the University of Q5 Antioquia.

Conflicts of interest
The authors declare that they have no conflicts of interest.