Investigation of SARS-CoV-2 lineages and mutations circulating in a university-affiliated hospital in South Korea analyzed using Oxford Nanopore MinION sequencing
Article information
Abstract
Objectives
Despite the introduction of vaccines, treatments, and massive diagnostic testing, the evolution of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has continued to overcome barriers that had slowed its previous spread. As the virus evolves towards increasing fitness, it is critical to continue monitoring the occurrence of new mutations that could evade human efforts to control them.
Methods
We performed whole-genome sequencing using Oxford Nanopore MinION sequencing on 58 SARS-CoV-2 isolates collected during the ongoing coronavirus disease 2019 pandemic at a tertiary hospital in South Korea and tracked the emergence of mutations responsible for massive spikes in South Korea.
Results
The differences among lineages were more pronounced in the spike gene, especially in the receptor-binding domain (RBD), than in other genes. Those RBD mutations could compromise neutralization by antibodies elicited by vaccination or previous infections. We also reported multiple incidences of Omicron variants carrying mutations that could impair the diagnostic sensitivity of reverse transcription-polymerase chain reaction-based testing.
Conclusion
These results provide an understanding of the temporal changes of variants and mutations that have been circulating in South Korea and their potential impacts on antigenicity, therapeutics, and diagnostic escape of the virus. We also showed that the utilization of the nanopore sequencing platform and the ARTIC workflow can provide convenient and accurate SARS-CoV-2 genomic surveillance even at a single hospital.
Introduction
Since its first emergence in December 2019, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has rapidly spread globally, and 4 months later the World Health Organization officially declared it a pandemic [1]. Since then, over 500 million human cases, including over 6 million deaths, have been reported worldwide so far, yielding several variants of interest (VOIs) and variants of concern (VOCs). The first case in South Korea was reported in January 2020, followed by a rapid increase in reported cases, which led to the first wave peak in February 2020 [2]. There have been several rises in the number of new infections in South Korea [2], including the latest surge of Omicron subvariants. Despite the strict restrictions on human activity and the introduction of vaccines and treatments, the evolution of the virus has continued to overcome barriers that had previously slowed its spread. The highest number of infections caused by new emerging variants has always exceeded the previous high, and the recent Omicron peak reached about 80 times higher than the previous Delta peak in South Korea [3].
During the course of the pandemic, more than 10 billion genome sequences of SARS-CoV-2 have been gathered worldwide, leading to the discovery that its continuous evolution is the fundamental driving force behind the prolonged pandemic [4]. Several genomic surveillance studies have also revealed that the accumulation of mutations has been heterogeneous along the genome and multiple genetic sites, such as the genes coding for the S and N proteins, were under positive selection [5,6]. As the virus evolves towards increasing fitness, it is necessary to continue monitoring the occurrence of new mutations that could increase infectivity or transmission efficiency, especially in response to these human efforts to control them.
In this study, we performed whole-genome sequencing using Oxford Nanopore MinION on 58 SARS-CoV-2 isolates collected during the ongoing coronavirus disease 2019 (COVID-19) pandemic (the 2nd, 4th, and 5th waves of infections) at a university-affiliated hospital in South Korea and tracked the emergence of mutations of certain variants responsible for the massive spikes in South Korea. Our results provide an understanding of the temporal change of variants and mutations that have been circulating in South Korea and the potential contribution of the evolving mutations on antigenicity, therapeutics, and diagnostic escape of the virus. We were also able to demonstrate that the utilization of the nanopore sequencing platform and the ARTIC workflow could provide convenient and accurate SARS-CoV-2 genomic surveillance even at a single hospital.
Materials and Methods
Sample Collection and Diagnostics
Nasopharyngeal and oral swab specimens from suspected cases were collected as part of routine surveillance at Kangdong Sacred Heart Hospital, Seoul, South Korea. RNA was extracted using a STARMag Universal Cartridge Kit (Seegene Inc., Seoul, Korea), and real-time reverse-transcription polymerase chain reaction (RT-PCR) with Seegene Allplex SARS-CoV-2 assay (Seegene Inc.) on a Bio-Rad CFX96 thermal cycler (Bio-Rad Laboratories, Hercules, CA, USA) was performed. Among samples that tested positive, those collected from June 2020 to September 2020 (Wuhan-type epidemic), November 2021 to December 2021 (Delta epidemic), and February 2022 to April 2022 (Omicron epidemic), which coincided with the 2nd, 4th, and 5th waves of the pandemic in South Korea, were subjected to whole-genome sequencing (Figure 1; Table S1).
Daily new confirmed cases and time series of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) waves in South Korea. Daily cases of SARS-CoV-2 in South Korea between January 20, 2020 and May 27, 2022 and 5 large waves are shown. The specimens used for analysis in this study were collected during the 2nd, 4th, and 5th waves, which are highlighted with asterisks.
Amplicon-Based Whole-Genome Sequencing
We performed whole-genome sequencing on 58 SARS-CoV-2 isolates following the ARTIC nCoV-2019 sequencing protocol. The cDNA synthesis and PCR amplification were performed according to the ARTIC amplicon sequencing protocol, using cDNA synthesis Automix (Seegene Inc.) in combination with the ARTIC nCoV-2019 V3 amplicon set. The PCR products were treated with a ligation kit (SQK-LSK109; Oxford Nanopore Technologies, Oxford, UK) and a barcoding kit (EXP-SFB001; Oxford Nanopore Technologies). Libraries were sequenced with the Oxford Nanopore MinION device using ONT R9.4.1 flow cells (Oxford Nanopore Technologies) for 12 hours.
Generation of Consensus Genomes and Mutational Profiling
For the construction of consensus sequences, the recommended ARTIC bioinformatics workflow (https://github.com/artic-network/fieldbioinformatics) was used. Generated sequences were manually assessed to verify the mutations called by the ARTIC pipeline. The genomes were classified into PANGO lineages using Pangolin v.4.0.6 [7], and mutation annotation was performed using Nextclade (https://clades.nextstrain.org) and Coronapp (http://giorgilab.unibo.it/coronannotator/). The genome sequences used for comparison were retrieved from the GISAID database (https://www.gisaid.org/) on April 6, 2022. Sequences were aligned using ClustalW 1.2.1 [8], and a maximum-likelihood phylogenetic tree was constructed using FastTree 2.1 [9] with 1,000 bootstrap replications. The Wuhan-Hu-1 sequence (RefSeq: NC_045512.2) was used as a reference genome and for rooting the tree as an outgroup. The tree was visualized using Interactive Tree of Life ver. 5 (https://itol.embl.de). Antibody escape was calculated using an escape calculator for the SARS-CoV-2 receptor-binding domain (RBD) [10].
Ethics Approval
This study was approved by the Institutional Review Board (IRB) of Kangdong Sacred Heart Hospital, Seoul, South Korea (IRB No: NON2021-001-001) and performed in accordance with the principles of the Declaration of Helsinki.
Results
Temporal Changes of SARS-CoV-2 Variants in South Korea
The ARTIC protocol generated 250,686±45,498 raw reads, totaling approximately 109±19 Mbp of sequences for each sample. The mean read qualities ranged between 11.6 and 13.9 (Table S1). All 58 sequenced samples showed at least 85% of the genome covered at least 20-fold, and 48 samples had more than 95% coverage (Table S2). Fifty-seven samples with at least 90% coverage were used in the subsequent whole-genome comparison (Table S1). In total, 6 SARS-CoV-2 lineages were identified in the 58 samples based on Pangolin [7] and Nextclade [11]. The samples collected between June 2020 and September 2020 (second wave of infections) were assigned to the B.1.497 (20C) lineage, and the samples collected between November 2021 and December 2021 (4th wave) belonged to the Delta VOC, which consisted of AY.69 (21I) and AY.122.5 (21J) (Table S1). Twenty-eight samples collected between March 2022 and April 2022 (5th wave) belonged to the Omicron VOC, including 11, 8, and 9 classified as BA.1.1 (21K), BA.2 (21L), and BA.2.3 (21L), respectively. The BA.2 and BA.2.3 lineages are currently the major causes of ongoing cases in South Korea [3]. The B.1.497 lineage had an average of 14.7 mutations, and the number of synonymous and non-synonymous mutations was similar on average (Figure S1; Table S3). However, an average of 45.1 and 44.4 mutations were identified in the 2 Delta subvariants, and the number of non-synonymous mutations was more than twice greater than that of synonymous mutations (Figure S1; Table S3). For the 3 Omicron subvariants (8 intact BA.1.1, 3 BA.2, and 5 BA.2.3 sequences were considered), the average numbers of mutations were 61.9, 66.7, and 70 for BA.1.1, BA.2, and BA.2.3 respectively, and the differences among them were also mainly due to the accumulation of non-synonymous mutations (42.1, 47, and 49.4 on average) (Figure S1; Table S3). The differences in the accumulation of non-synonymous and synonymous mutations were greater in the S region, compared to the other gene regions (Figure S1). As the virus evolved, while the number of synonymous mutations maintained below 7 in the S gene, non-synonymous mutations markedly increased to an average of 28, 22, and 23 in BA.1.1, BA.2, and BA.2.3, respectively (Figure S1).
The divergence between lineages was the greatest in the spike protein. The BA.1.1 subvariant contained about 32 new mutations in the spike encoding region compared to the Delta variants, and the BA.2 and BA.2.3 subvariants gained about 9 new mutations (T19I, A27S, G142D, S371F, V213G, T376A, D405N, R408S, and 24–26 deletion) compared to BA.1.1 (Tables S2, S3). Meanwhile, the 2 subvariants lacked 14 mutations, including 69–70 del, 142–144 del, and 211 del. Although most of the amino acid changes in the spike were fixed at the individual level, KD-Cov-O19 carried 2 unique mutations, L5F and H49Y (Tables S2, S3).
Impacts of the Emergence of Mutations on Infectivity and Therapeutics
As evident from Figure 2 and Table S4, the differences among lineages were more pronounced in the spike gene. While the B.1.497 possessed only 1 fixed mutation in the gene, 6–7 and 27–31 fixed mutations were identified in genomes of the Delta (AY.122.5 and AY.69) and the Omicron (BA.1.1, BA.2, and BA.2.3) variants, respectively. Markedly, the majority of the mutations gained by the Omicron variants were located in the RBD region of the gene, which had 7.5–8 times more substitutions (16 for BA.1.1 and 15 for BA.2 and BA.2.3) than the Delta subvariants in the region, and all these substitutions were non-synonymous changes. Of the 16 changes in Omicron BA.1.1 and 15 changes in the Omicron BA.2 and BA.2.3 subvariants in the RBD, 8 (K417N, G446S, E484A, Q493R, G496S, Q498R, N501Y, and Y505H) and 6 (K417N, E484A, Q493R, Q498R, N501Y, and Y505H), respectively, mapped to the positions proposed to be critical for binding to human ACE2 [12], and the residues remained unchanged in other lineages (Figure 3A).
Genome-wide mutational profiling of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) isolates. This figure illustrates genome-wide patterns of mutations in 6 SARS-CoV-2 lineages responsible for the major waves in South Korea. A maximum-likelihood phylogenetic tree was constructed to show the relationship between the isolates, and the tree labels were colored based on the lineages. Each colored circle and square represent different types of mutations, located according to their actual position in each genome.
Mutations in the spike protein receptor-binding domain (RBD) may affect binding to human ACE2 and antibodies. (A) Differences in amino acid residues of known binding sites to human ACE2. (B–E) Antibody escape map for 33 neutralizing antibodies targeting RBD. The blue line shows escape mediated by mutations at each site, and the gray line represents the original escape map in the absence of mutations. The orange circles represent each mutation. The bar graphs show the total antibody binding remaining after combinations of mutations occurred. (B) B.1.497, (C) AY.122.5 and AY.69, (D) BA.1.1, (E) BA.2 and BA.2.3.
NTD, N-terminal domain; FP, fusion peptide; HR1, heptad repeat 1; HR2, heptad repeat 2; TM, transmembrane domain.
The antibody-escape calculator [10] estimated that the combination of the RBD mutations of the Delta subvariants reduced the fraction of binding antibodies to below 0.9 (Figure 3C). The mutations that emerged in the Omicron variants could lead to more immune escape by dropping the fraction to 0.5 (BA.1.1) and 0.6 (BA.2 and BA.2.3), and the difference between these subvariants was induced by additional R346K, N440K, G446S, and G496S mutations in the BA.1.1 lineage (Figure 3D, E).
In the Mpro coding region, which is an important antiviral drug target of Paxlovid (Pfizer) for treatment of the virus, 1 B.1.497 isolate (KD-Cov-N13) was observed to have a change in the 168 position, and 2 AY.69 isolates, D7 and D11, harbored M264I and a synonymous change in the 10078 nucleotide position, respectively (Table S3). The BA.1.1 lineage carried 1 mutation (P132H), and 4 types of substitution (synonymous changes in the 10198, 10432, and 10447 nucleotide positions and P132H) were observed in BA.2 and BA.2.3. The substitution at the 10432 position occurred only in KD-Cov-O11, BA.2.3, which had the most substitutions in the nsp5 region among our isolates.
The Effect of the Emergence of Mutations on Diagnostic Sensitivity
Routine diagnostic testing has been performed at our diagnostic facility since the beginning of the SARS-CoV-2 pandemic. However, since February 2022, we have observed a few cases (KD-CoV-25, 26, 27, 28, and 29) showing discrepant results in target gene detection with a commercial assay (Allplex SARS-CoV-2 assay) which targets N, E, and RNA-dependent RNA polymerase (RdRP)/S gene fragments. KD-Cov-O25 and O26, a mother-child pair, showed N gene drop-out, despite having cycle threshold (Ct) values below 25 for both the E and RdRP/S genes (Table 1). In the KD-Cov-O27 and O28 samples, another mother-child pair, the Ct values for N and the other 2 (E and RdRP/S) differed by about 5 cycles, and regarding KD-Cov-O29, the Ct value for RdRP/S was 8 cycles higher than those of the E and N genes (Table 1). Deletions were identified in position 28877-28894, leading to the deletion of 6 amino acids (203–208 del) in the N protein of the KD-CoV-O25 and O26 pair. The KD-CoV-O27 and O28 pair carried an A208V mutation (C28896T) in the protein, overlapping the deleted region of the above-mentioned isolates. The only change observed in the KD-Cov-O29 genome was an additional M657T mutation (T15437C) in NSP12 compared to BA.2.3 isolates with normal Ct values. Although information on the primers and probes used in commercial SARS-CoV-2 assays is not generally disclosed, the mutations observed in the above cases may affect the RT-PCR diagnostic efficacy of the Allplex SARS-CoV-2 assay. To date, a total of 1,281 cases with deletions in the 203-208 position were reported on GISAID (10,140,904 deposited N gene sequences as of April 6, 2022). Interestingly, this deletion has been sporadically occurring in Canada, the USA, Brazil, Germany, Italy, and several other countries since its first appearing on GISAID in early 2020 (EPI_ISL_3717678, Isolation: India, March 19, 2020) before the emergence of Omicron. Furthermore, it was not specific to the Omicron variant and was also observed in Alpha, Delta, Gamma, Omicron, and other lineages. It was most prevalent in the USA and Brazil (397 and 186 instances, respectively). However, in South Korea, the deletion was first observed in January 2022 in the genome sequence of the Omicron BA.1.1 variant on GISAID and has maintained a low frequency (7 instances/11,785 Omicron sequences). It has also been exclusively observed in the Omicron subvariant (BA.1.1) so far in South Korea. The A208V mutation has been observed in 4,447 genome sequences, and the majority of them originated from the USA (1,009) and England (897).
Ct values for specimens with the Allplex SARS-CoV-2 assay and potentially associated mutations
Discussion
In this study, we performed whole-genome sequencing on 58 SARS-CoV-2 isolates collected during the ongoing COVID-19 pandemic at a single hospital to track the emergence of mutations of certain variants responsible for the massive spikes in South Korea. The differences among lineages were more pronounced in the spike gene, especially in the RBD region. Since changes in the RBD region may substantially change the antigenicity and susceptibility to pre-existing antibodies [12,13], mutations of emerging viruses in this region should be closely monitored. The new mutations gained by the Omicron variants were mostly located in the critical binding sites to human ACE2. Interestingly, G446S and G496S were exclusively detected in BA.1.1. Another study reported that G496 was the residue introducing the most significant effects on the RBD-ACE2 interaction [14], and the G496S mutation was also suggested to reinforce or generate chemical bonding at the interface, along with G446S and N501Y [15]. Of the mutually detected mutations in the 3 Omicron subvariants, a previous study showed that Q493R and N501Y, along with S477N, formed a new interaction with human ACE2, which enhanced the binding of the Omicron variant to human ACE2 [16]. The Y505H mutation did not seem to have an impact on binding [17].
We then analyzed how the emerging RBD mutations could compromise neutralization by antibodies elicited by vaccination or natural infection, using a recently developed escape calculator tool that utilized a deep mutational scanning approach. Previous studies demonstrated that the calculator correlated well with experimentally measured neutralization titers [10]. The results indicated that the combination of RBD mutations that emerged in the Omicron variants could lead to more immune escape compared to the previous lineages. The difference between the Omicron subvariants was induced by additional R346K, N440K, G446S, and G496S mutations in the BA.1.1 lineage (Figure 3D, E). Motozono et al. [18] showed that Omicron BA.1 with the G446S mutation, located adjacent to NF9 epitope, showed greater inhibition of replication by vaccine-induced T cells as a result of enhanced presentation and recognition of the antigenic peptide on the spike protein. As shown in Figure 3B–E, the mutation at position 484 was among the largest drivers of antigenic changes [19]. In a study using 19 monoclonal antibodies, substitutions at this position occurred the most frequently [20], and all amino acid changes to A, K, Q, P, D, and G reduced recognition by antibodies [20,21]. The E484A mutation detected in Omicron variants also ablated neutralization by bamlanivimab, a therapeutic monoclonal antibody [22]. The E484K mutation carried in Beta, Gamma, Zeta, Eta, and Theta variants [23] is also well known for its capability of reducing antibody binding, and it emerged in a patient during the use of monoclonal antibodies or convalescent plasma therapy [24].
Rapid emergence of new variants with higher mutation loads also raises concern for the clinical effectiveness of currently used treatments. The non-monoclonal antibody therapeutics currently approved and used to treat SARS-CoV-2 are nirmatrelvir (an antiviral component of Paxlovid), remdesivir, and molnupiravir. Nirmatrelvir inhibits Mpro, and both remdesivir and molnupiravir are RdRP inhibitors. According to previous studies, there are 26 binding site amino acid residues in NSP5 with nirmatrelvir, and key interaction sites were discovered [25,26]. Interestingly, the P168S mutation identified only in KD-Cov-N13 was located at this site, but the sample was collected in September 2020, which indicates that the mutation was not driven by selective drug pressure. Although other substitutions were not located in the direct binding site, the synonymous change in the 10198 nucleotide position, frequently detected in lineages BA.2 and BA.2.3, was located in the proximity of the S2 binding pocket of nirmatrelvir [25]. P132H, one of the most prevalent NSP5 variants, was also detected in all the Omicron isolates in this study, but several studies have shown that it is unlikely to affect nirmatrelvir binding affinity, as the residue is located in a flexible turn [27,28]. We speculate that no mutations affecting Paxlovid resistance have occurred in this study period. In addition, no mutations were detected in RdRP residues that interact with remdesivir or molnupiravir [29,30].
Paxlovid has been globally used since it was authorized by the Food and Drug Administration in December 2021. Research has shown that although certain mutations such as P132H, G15S, and K90R have become frequently reported, the binding and inhibitory effects of nirmatrelvir were still not compromised by those variants [31]. Nonetheless, it is obvious that, with the advent of the Omicron variant, the frequency and number of substitutions in the NSP5 region have markedly increased compared to earlier VOCs and VOIs, and the fact that several mutational changes have been occurring at the inhibitor binding sites, albeit at a relatively low frequency [26], indicates that continuous genomic monitoring and characterization of those mutations in the protease sequence should be performed.
After the Omicron variant emerged in South Korea, we observed a few mutations affecting the detection of the N gene. By carrying these mutations, the virus may have the ability to evade detection, thus having selective advantages such as increased transmission [32]. Recurrence of the mutations in the same genetic region across various lineages during the course of the SARS-CoV-2 pandemic also implies certain beneficial roles of the mutation. The 203−208 del mutation is positioned in the linker region, which links the N-terminal and C-terminal of the N protein. Previous studies have suggested that this deletion might impact phase separation, promoting the function of the region and possibly provides advantages in the packaging and replication process of the virus [33−35]. The N gene was generally considered a highly conserved region and has been used for the target region of several RT-PCR diagnostic kits [36]. However, recent genomic surveillance studies have shown that this gene, particularly the linker region, is one of the least conserved regions in the SARS-CoV-2 genome [33,35,37].
There also have been reports of similar detection failure problems in several diagnostic kits. A mutation in the E gene (C26340T) was reported to affect E gene detection by the Cobas SARS-CoV-2 assay (Roche Inc., Basel, Switzerland) [38] in 2020. Furthermore, C29200A [39] and G29179T [40] mutations in the N gene were observed to impair the sensitivity of Xpert Xpress SARS-CoV-2 assay (Cepheid Inc., Sunnyvale, CA, USA) in 2020. In early 2021, Zannoli et al. [32] collected VOC B.1.17 specimens showing N gene drop-out using the Allplex SARS-CoV-2 assay, which resulted from the 207−208 deletion or 208-209 indel. Furthermore, more than 100 specimens, mostly Gamma variants, with impaired N gene detection by the GeneFinder Kit were reported by Lesbon et al. [33], and 3 different mutations in the gene were identified as associated mutations, including 202−207 deletion, between April 2020 and July 2021. The constant occurrence of these mutations could affect the diagnostic sensitivity of RT-PCR-based testing and lead to an increase in false-negative results, especially in samples with a low viral load or in single target assays. Therefore, continuous monitoring of the mutations and the consequent sensitivity of various diagnostic kits is crucial for maintaining accurate diagnostics.
The results provided in this study will contribute to an understanding of the potential impacts of the emerging mutations on evading global efforts to control the virus, thus underlining the need for continual genomic surveillance and characterization of the evolution of SARS-CoV-2. Although whole-genome sequencing plays a key role in pandemic surveillance, the clinical adoption of sequencing has been limited by several aspects, such as higher costs, technical and bioinformatics expertise, and time for implementation. In this respect, the present study showed that utilization of the nanopore sequencing platform and the ARTIC workflow could offer rapid, cost-effective, and reliable sequencing of SARS-CoV-2 samples even at a single hospital. This strategy can serve as a method for increasing the national genomic surveillance capacity for SARS-CoV-2 and other future pathogens in pandemics and epidemics by enabling private labs, academia, and other regional institutions to actively engage in surveillance.
Supplementary Material
Supplementary data are available at https://doi.org/10.24171/j.phrp.2022.0183.
Information on the SARS-CoV-2 specimens used in this study and sequencing statistics;
Coverage depths of each amplicon;
Mutations that occurred in each genome of 58 SARS-CoV-2 isolates;
Distribution of mutations in the spike proteins of 58 SARS-CoV-2 isolates. Positions that are highlighted in yellow represent regions with low sequencing depths.
Comparison of the number of non-synonymous and synonymous mutations in each genetic region.
Notes
Ethics Approval
This study was approved by the Institutional Review Board (IRB) of Kangdong Sacred Heart Hospital, Seoul, South Korea (IRB No: NON2021-001-001) and performed in accordance with the principles of the Declaration of Helsinki.
Conflicts of Interest
The authors have no conflicts of interest to declare.
Funding
This work was supported by Kangdong Sacred Heart Hospital Fund (2021-05) and the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI; grant number: HI20C0071), funded by the Ministry of Health & Welfare, Republic of Korea.
Availability of Data
The datasets generated during the current study are available in the NCBI GenBank accession no. ON797335-ON797392.
Authors’ Contributions
Conceptualization: JSK, HK; Data curation: HK; Formal analysis: HK, HySK; Funding acquisition: JSK; Investigation: SHC, HySK; Methodology: SHC, JSK; Project administration: JSK; Resources: JSK; Software: HK; Supervision: KHH; Validation: HaSK, WS; Visualization: HK; Writing–original draft: HK, JSK, SHC, HySK, HaSK, WS, KHH; Writing–review & editing: all authors.