Genetic Diversity of SARS-CoV-2 over a One-Year Period of the COVID-19 Pandemic: A Global Perspective

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) caused a global pandemic of coronavirus disease in 2019 (COVID-19). Genome surveillance is a key method to track the spread of SARS-CoV-2 variants. Genetic diversity and evolution of SARS-CoV-2 were analyzed based on 260,673 whole-genome sequences, which were sampled from 62 countries between 24 December 2019 and 12 January 2021. We found that amino acid (AA) substitutions were observed in all SARS-CoV-2 proteins, and the top six proteins with the highest substitution rates were ORF10, nucleocapsid, ORF3a, spike glycoprotein, RNA-dependent RNA polymerase, and ORF8. Among 25,629 amino acid substitutions at 8484 polymorphic sites across the coding region of the SARS-CoV-2 genome, the D614G (93.88%) variant in spike and the P323L (93.74%) variant in RNA-dependent RNA polymerase were the dominant variants on six continents. As of January 2021, the genomic sequences of SARS-CoV-2 could be divided into at least 12 different clades. Distributions of SARS-CoV-2 clades were featured with temporal and geographical dynamics on six continents. Overall, this large-scale analysis provides a detailed mapping of SARS-CoV-2 variants in different geographic areas at different time points, highlighting the importance of evaluating highly prevalent variants in the development of SARS-CoV-2 antiviral drugs and vaccines.


Introduction
Since December 2019, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused a global pandemic. As of 16 January 2021, more than 93 million cases and 2 million deaths have been reported, according to the daily update of the World Health Organization. As a virus of the Coronaviridae family, SARS-CoV-2 is a positive-sense, singlestranded RNA virus whose genome contains approximately 30,000 nucleotides, making it one of the largest genomes among RNA viruses. The coding region of the SARS-CoV-2 genome encodes for nonstructural proteins (NSPs) and structural proteins in different open reading frames (ORFs). The genetic makeup of polyproteins pp1a and pp1ab, encoded by ORFs 1a and 1b at the 5 -end, contains 16 nonstructural proteins (numbered from NSP1 to NSP16). The 3 -end of the genome (21 to 29 kb) encodes for at least 6 accessory proteins (ORF3a, ORF6, ORF7a, ORF7b, ORF8, ORF10) [1] and 4 structural proteins: spike glycoprotein (S), envelope protein (E), membrane glycoprotein (M), and nucleocapsid protein (N) ( Figure 1A). Since viral proteins such as RNA-dependent RNA polymerase and spike glycoprotein are key drug targets [2], genetic variations in the genome may exert an important impact on the efficacy of current vaccines and antiviral treatments [3]. Global surveillance of SARS-CoV-2 genetic variants is important to control the coronavirus disease 2019 (COVID- 19) pandemic. Previous analyses of SARS-CoV-2 genomes have characterized proximal origin [4][5][6], evolutionary rates [7], spatiotemporal distributions [8], transmission dynamics in many countries [9][10][11][12][13][14][15], phylogenetic networks of dominant strains [16][17][18][19], intragenomic diversity [20], epistatic interactions between viral genes [21], genetic variants in response to disease severity and immune responses [22][23][24][25], and the impact of spike variants on current vaccines and antibodies [26][27][28][29][30]. Despite its long genome structure, SARS-CoV-2 encodes a 3′-to-5′ exoribonuclease called NSP14 that provides the proofreading activity to correct errors made by the viral RNA-dependent RNA polymerase, thereby promoting a high fidelity of genome replication [31,32]. Nevertheless, some variants, such as D614G in spike, can escape proofreading and significantly shape the global dispersal of SARS-CoV-2 infections [33]. Moreover, circulating variants may exert a potential impact on the development of vaccines and antiviral drugs to control the COVID-19 pandemic [34]. Due to the dynamic dispersal of viral strains worldwide, it remains important to continue the surveillance of SARS-CoV-2 variants from a global perspective [35].
Here, we have conducted global surveillance of circulating variants in SARS-CoV-2 based on a large-scale dataset of whole-genome sequences from 62 countries worldwide. For the development of novel vaccines and antiviral drugs, this study contributes to a deeper understanding of diversity and global surveillance of key viral variants of SARS-CoV-2. Global surveillance of SARS-CoV-2 genetic variants is important to control the coronavirus disease 2019 (COVID-19) pandemic. Previous analyses of SARS-CoV-2 genomes have characterized proximal origin [4][5][6], evolutionary rates [7], spatiotemporal distributions [8], transmission dynamics in many countries [9][10][11][12][13][14][15], phylogenetic networks of dominant strains [16][17][18][19], intragenomic diversity [20], epistatic interactions between viral genes [21], genetic variants in response to disease severity and immune responses [22][23][24][25], and the impact of spike variants on current vaccines and antibodies [26][27][28][29][30]. Despite its long genome structure, SARS-CoV-2 encodes a 3 -to-5 exoribonuclease called NSP14 that provides the proofreading activity to correct errors made by the viral RNA-dependent RNA polymerase, thereby promoting a high fidelity of genome replication [31,32]. Nevertheless, some variants, such as D614G in spike, can escape proofreading and significantly shape the global dispersal of SARS-CoV-2 infections [33]. Moreover, circulating variants may exert a potential impact on the development of vaccines and antiviral drugs to control the COVID-19 pandemic [34]. Due to the dynamic dispersal of viral strains worldwide, it remains important to continue the surveillance of SARS-CoV-2 variants from a global perspective [35].
Here, we have conducted global surveillance of circulating variants in SARS-CoV-2 based on a large-scale dataset of whole-genome sequences from 62 countries worldwide. For the development of novel vaccines and antiviral drugs, this study contributes to a deeper understanding of diversity and global surveillance of key viral variants of SARS-CoV-2.

Data Collection
On 16 January 2021, we extracted SARS-CoV-2 complete genome sequences with high coverage (with <1% undefined bases and <0.05% unique amino acid (AA) substitutions and no insertion/deletion unless verified by submitter) from the GISAID database (https: //www.gisaid.org/, the accessed date: 16 January 2021). We also downloaded the sequence metadata regarding the geographical locations and collection times. We acknowledge all contributors who have kindly deposited and shared their genome data on GISAID. A list of genome sequence acknowledgments is provided in Data S1. The complete genome of the Wuhan-Hu-1 isolate (NCBI accession NC_045512) was also collected as a reference strain.
Sequence quality control was conducted as follows: (i) one sequence per patient was sampled, while multiple sequences sampled from one patient were screened using the metadata regarding patient ID, authors, lab origin, and country origin; and (ii) sequences from nonhuman origins were removed.

Multiple Sequence Alignment
Nucleotide sequences were aligned to the reference sequence (Wuhan-Hu-1) using the option of "addfragments" in MAFFT version 7.471 [36], and output alignments were manually scrutinized by Seaview version 5.0 (http://doua.prabi.fr/software/seaview, the accessed date: 16 January 2021). The 5 and 3 untranslated regions were trimmed, and only sequences that contained full-length coding regions were retained. Inhouse MATLAB code (available at https://github.com/MiaoMiaorrk/sequence_analysis, the accessed date: 16 January 2021) was prepared to extract the nucleotide regions of 26 viral proteins, while the coding region of the viral polymerase NSP12 was concatenated by considering programmed −1 ribosomal frameshifting. The corresponding amino acid sequences were subsequently obtained for analyzing protein diversity and amino acid variants.

Protein Sequence Diversity
Polymorphic sites are amino acid or nucleotide positions that have at least one substitution compared with the reference Wuhan-Hu-1, ignoring gaps and ambiguities [30]. Given multiple sequence alignments, the amino acid substitution rate was defined as: Setequal(a, b) = 1 , i f a = b and a, b are both valid; 0 , otherwise.
A(i, j) is the element at the i position of the j th sequence, and A(i, 0) is the element at the i position of the reference Wuhan-Hu-1. L is the length of the protein, and N is the number of sequences except for the reference sequence. N Invalid is the total number of invalid positions containing gaps or ambiguous letters. Setequal(A(i, j), A(i, 0)) is the function that checks if two elements A(i, j) and A(i, 0) are equal.

Selective Pressure Analysis
To analyze the selective pressure on the protein-coding regions of SARS-CoV-2, the dN/dS ratio that measures the nonsynonymous (dN) to synonymous (dS) substitution rates at the coding region was calculated using KaKs_Calculator 2.0 with the method of modified YN (MYN) and the genetic code in Table 1 (standard code) [37]. The positive selection was defined by dN/dS > 1.

Phylogenetic Analysis
At each sample collection date, from 24 December 2019 to 12 January 2021, three genomic sequences were randomly subsampled for each continent, while the sequences were all collected if the number of sequences was less than 3. This led to a subset of SARS-CoV-2 genome sequences (n = 5488). Multiple sequence alignments of this subset were subsequently imported to reconstruct a maximum-likelihood phylogenetic tree using NextStrain tools (https://nextstrain.org, the accessed date: 16 January 2021).

Basic Characteristics of SARS-CoV-2 Genome Sequences
Our study analyzed a large-scale dataset of SARS-CoV-2 genome sequences (n = 260,673) that fulfilled the quality control criteria. During the pandemic of COVID-19, from 24 December 2019 to 12 January 2021, these full-length sequences were sampled from 62 countries located in Europe (n = 166,443), North America (n = 61,232), Asia (n = 13,709), Oceania (n = 13,166), South America (n = 3216), and Africa (n = 2907). The collection time of SARS-CoV-2 genome sequences is visualized in Figure 1B. Most sequences were collected during the first wave (March to May 2020) before the summer in the Northern Hemisphere and during the second wave (September 2020 to January 2021) in Europe and North America, where the COVID-19 pandemic remains severe.
Due to the worldwide spread of SARS-CoV-2 infections, dynamic patterns of protein diversity were observed over time, especially in six viral proteins: ORF10, nucleocapsid, ORF3a, spike, NSP12, and ORF8 ( Figure 2A). Unlike other proteins, whose protein diversities increase gradually over time, the protein diversity of ORF10 has increased sharply since July 2020. Figure 2B shows that the highest substitution rate of the ORF10 protein appeared in Europe, followed by Africa. Although the increasing patterns of protein diversity of the remaining four proteins, from December 2019 to January 2021, were roughly the same, the temporal patterns of spike and ORF10 diversity were different between Africa and other continents ( Figure 2B). Protein diversity of spike and ORF10 reached their peaks in November and December 2020, respectively.

Positive Selection is Rare in SARS-CoV-2 Protein-Coding Genes
To analyze selective pressure on SARS-CoV-2 protein-coding regions, we estimated the dN/dS ratios of nonsynonymous and synonymous substitution rates ( Figure S1). The median dN/dS ratio for all protein-coding genes was 0.31 (IQR: 0.24 to 0.65). All viral protein-coding genes had a median dN/dS ratio less than 1, and the highest level of median dN/dS ratio was observed in ORF9 (0.65, IQR: 0.49 to 0.65), followed by ORF6 (0.56) and ORF7b (0.50) (Table S1).

Positive Selection Is Rare in SARS-CoV-2 Protein-Coding Genes
To analyze selective pressure on SARS-CoV-2 protein-coding regions, we estimated the dN/dS ratios of nonsynonymous and synonymous substitution rates ( Figure S1). The median dN/dS ratio for all protein-coding genes was 0.31 (IQR: 0.24 to 0.65). All viral protein-coding genes had a median dN/dS ratio less than 1, and the highest level of median dN/dS ratio was observed in ORF9 (0.65, IQR: 0.49 to 0.65), followed by ORF6 (0.56) and ORF7b (0.50) (Table S1).
We next analyzed the amino acid variants of SARS-CoV-2. There were 8484 (86.95%) polymorphic sites calculated as the number of sites having one or more substitutions compared to the reference sequence (Wuhan-Hu-1, NCBI accession NC_045512) across the 9757 amino acid positions in 26 concatenated proteins. In total, 25,629 substitutions were identified at 8484 amino acid positions. Notably, 587 of 1273 fully conserved positions were mostly located in the enzymatic proteins of NSP12 (17.05%), NSP3 (15.16%), and spike (13.90%). The frequency of variants on different continents was calculated by dividing the number of sequences carrying a given variant on a specific continent by the total number of sequences on that continent; 115 polymorphic sites occurred with a frequency greater than 0.5%. Highly polymorphic sites that occurred with a frequency greater than 1% were mostly observed in nucleocapsid (3.10%), ORF10 (2.63%), ORF8 (2.48%), ORF3a (2.18%), spike (0.71%), and NSP12 (0.64%).
We next analyzed the amino acid variants of SARS-CoV-2. There were 8484 (86.95%) polymorphic sites calculated as the number of sites having one or more substitutions compared to the reference sequence (Wuhan-Hu-1, NCBI accession NC_045512) across the 9757 amino acid positions in 26 concatenated proteins. In total, 25,629 substitutions were identified at 8484 amino acid positions. Notably, 587 of 1273 fully conserved positions were mostly located in the enzymatic proteins of NSP12 (17.05%), NSP3 (15.16%), and spike (13.90%). The frequency of variants on different continents was calculated by dividing the number of sequences carrying a given variant on a specific continent by the total number of sequences on that continent; 115 polymorphic sites occurred with a frequency greater than 0.5%. Highly polymorphic sites that occurred with a frequency greater than 1% were mostly observed in nucleocapsid (3.10%), ORF10 (2.63%), ORF8 (2.48%), ORF3a (2.18%), spike (0.71%), and NSP12 (0.64%).
In addition, the major variants of each SARS-CoV-2 protein are mapped in Figure 3. For each position, the Wuhan-Hu-1 index is shown at the top, followed by the variants with a frequency >1%. Fifteen variants with the prevalence >5% were identified, including (i) four variants in nucleocapsid: S194L (6.29%), R203K (28.  For each site, the reference index is shown at the top, followed by variants with a frequency >1%. Variants highlighted with green superscripts were those with frequencies >5%. For each site, the reference index is shown at the top, followed by variants with a frequency >1%. Variants highlighted with green superscripts were those with frequencies >5%.
( Figure 5C), and V30L in ORF10 increased continuously in Europe and Africa, while their emergence in other continents was less common. (iii) Two neighboring variants in nucleocapsid, R203K and G204R ( Figure 5C), showed roughly similar frequencies over time in the same geographic areas. (iv) Temporal trends of the NSP2 variant T85I ( Figure 5D) and the ORF3a variant Q57H were distinct on different continents. From March to December 2020, the prevalence of T85I and Q57H was higher in North America than on other continents.

Discussion
This study presents global surveillance of SARS-CoV-2 genetic diversity and circulating variants based on a large-scale dataset of 260,673 sequences sampled from December 2019 to January 2021. Based on this first one-year survey of genetic diversity, our study describes the temporal and geographical trends of circulating variants on different continents. Despite the proofreading activity of viral exoribonuclease, an increasing prevalence of key genetic variants, such as D614G in viral spike and P323L in NSP12, was observed in many countries and continents. Early surveillance of potentially evolving viral strains with higher transmission fitness and greater infectivity plays an important role in controlling the ongoing COVID-19 pandemic.
Among 26 viral proteins, the highest protein sequence diversity was observed in spike, nucleocapsid, ORF3a, and ORF10 ( Table 1). As a critical target of many vaccines and neutralizing antibodies, the spike protein interacts with cellular receptors, such as angiotensin-converting enzyme 2 (ACE2), for viral entry [38,39]. Although whether the spike variants may weaken the performance of vaccines is still under evaluation [40,41], antibody therapeutics and vaccine development still need to consider the impact of spike variants on the antigenicity of the virus [42]. As an indispensable protein, nucleocapsid is the major target of primers and probes in many real-time reverse transcription polymerase chain reaction (RT-PCR) tests, while nucleotide substitutions in the primer binding region may affect the sensitivity of some diagnostic assays [43]. Although single nucleotide substitutions may not be an issue for solid assays that target multiple loci (n > 2) in the genome [44], close monitoring of SARS-CoV-2 nucleotide substitutions that are located within the primer and probe regions of RT-PCR assays remains important. ORF3a plays a role in virulence, ion channel formation, viral release, and apoptosis [45], while ORF10 is an accessory protein that is dispensable for viral replication in humans [46]. Our study presents that the highest substitution rate of the ORF10 protein appeared in Europe, followed by Africa ( Figure 2B). Previous studies have reported the import of SARS-CoV-2 from Europe to Africa [47,48], which explains some similar patterns between Europe and Africa. To support surveillance, a paucity of SARS-CoV-2 sequence data from Africa requires further sequencing efforts on the African continent [47]. Additionally, other viral proteins are also known to interact with many host proteins that participate in multiple biological processes such as protein trafficking, translation, and transcription [49]. Genetic variants that promote viral fitness and escape may exert a potential impact on current diagnostic tests, vaccines, antiviral strategies, and immune responses [23,40,50]. Close monitoring of key viral variants is critical for the development and optimization of vaccines and antiviral therapies, according to previous experience from the management of viral infections such as influenza virus, respiratory syncytial virus, and human immunodeficiency virus [51][52][53][54][55][56][57].
Highly prevalent variants have been observed in SARS-CoV-2. A high prevalence of important variants, such as D614G (93.88%) in spike and P323L (93.74%) in NSP12, was confirmed by our large-scale study. Compared with wildtype D614, the D614G variant changes the dynamic structures of the viral spike to bind with human angiotensinconverting enzyme 2 [58] and increases transmission fitness with enhanced viral loads in the upper respiratory tract of SARS-CoV-2 cases [34]. The D614G variant is located at a region interfacing with the neighboring subunit, and it potentially alters the trimer stability of the viral spike to enhance membrane fusion and host entry [59]. In addition, D614G increases the replication fitness of SARS-CoV-2 in cell cultures and enhances viral transmission in hamsters [60]. Furthermore, D614G in spike is coevolving with P323L in NSP12, with strong allelic associations [16]. The coexistence of D614G (spike), P323L (NSP12), and C241U (5 -UTR) may contribute to increased transmission fitness [61]. The P323L variant might alter the secondary and tertiary structures of NSP12 to interact with NSP8, thereby affecting viral replication [62,63]. Additionally, other variants like S25L in NSP7, S194L, R203K, G204R, and T205I in nucleocapsid, T85I and I120F in NSP2, S477N in spike, and Q57H in ORF3a have also been confirmed by previous studies [62,64]. Although many SARS-CoV-2 variants have been observed, their associations with viral loads and infectivity need more investigation. Vaccines may need to be updated periodically to avoid the potential impact of newly arising SARS-CoV-2 variants on the clinical efficacy of the vaccines [65].
Key SARS-CoV-2 variants in certain clades may play an important role in viral transmission and adaptability. Monitoring the dynamics of SARS-CoV-2 clades is crucial because different viral clades, such as B.1.1.7, may be associated with viral loads and disease severity [66]. Due to the evolving nature of SARS-CoV-2, the emergence of some variants in B.1.1.7, B.1.351, and P.1 clades may affect viral transmissibility and antigenic profiles [65,67]. Based on large-scale full-length genome sequences, our findings showed distinct geographical and temporal patterns of SARS-CoV-2 clades on different continents (Figure 7). We observed the dynamic patterns of SARS-CoV-2 clades over time; for instance, Clade 20E (EU1) has become the dominant group worldwide since mid-September 2020. The 20H/501Y.V2 and 20I/501Y.V1 clades were observed in Africa in October 2020 and in Europe in November 2020, respectively. SARS-CoV-2 variants in different clades might exert an impact on some vaccines and antibodies [27,65,68]. Spike variants such as D614G and N501Y, contained in the 20H/501Y.V2 and 20I/501Y.V1 clades, may lead to changes in viral antigenicity that are harmful to monoclonal antibody therapies and vaccine protection [67], thereby mediating potential escape from vaccine response [41]. Since a cure against COVID-19 is still lacking, the existence of current and novel circulating clades that are globally overdispersed highlights the importance of rapid interventions to reduce viral transmissions [10].
There are limitations in our study. First, our analyses focused on genetic diversity over the past one-year period, while global surveillance of new clades or lineages is still needed during the ongoing COVID-19 pandemic. Second, patient information on disease severity, immune responses, and transmission history was largely lacking in our retrieved dataset. Future studies are needed to reveal potential associations of genetic variants with disease progression and transmission fitness. Third, a limited number of full-length genome sequences from Africa and Asia were deposited in the public databases, which limited our analyses to fully characterize the genetic diversity of SARS-CoV-2 in Africa and Asia. Future analyses should integrate more sequences from different countries and continents.

Conclusions
Despite its global spread, SARS-CoV-2 harbors only a small proportion of highly prevalent variants that determine the classification of viral clades or lineages. Based on a large-scale dataset of 260,673 genomic sequences, our study presents a comprehensive mapping of highly prevalent SARS-CoV-2 variants as well as the dynamic changes of clades during the one-year pandemic of COVID-19. Genomic surveillance is important for monitoring the genetic variants that may be resistant to current antiviral drugs and vaccines [35,69]. To effectively control the COVID-19 pandemic, further studies are still needed to evaluate the impact of highly prevalent variants on antiviral drugs and vaccines before their wide application.