Genetic diversity of SARS-CoV-2 in South America: demographic history and structuration signals.

In 2020, the emergence of SARS-CoV-2 caused a global public health crisis with significant mortality rates and a large socioeconomic burden. The rapid spread of this new virus has led to the appearance of new variants, making the characterization and monitoring of genetic diversity necessary to understand the population dynamics and evolution of the virus. Here, a population-genetics-based study was performed starting with South American genome sequences available in the GISAID database to investigate the genetic diversity of SARS-CoV-2 on this continent and the evolutionary mechanisms that modulate it.

Introduction

Many human infectious diseases have emerged because of the pathogen’s ability to jump hosts [1, 2]; this is the case with betacoronaviruses [3]. Like SARS-CoV and MERS-CoV [4, 5], a new member of the genus Betacoronavirus termed SARS-CoV-2 [6] emerged as a zoonotic virus [7, 8]. It was recorded for the first time at the end of 2019 in Wuhan, China [9], and it is responsible for the COVID-19 pandemic. The genome of SARS-CoV-2 is a 29.9-kb, positive-sense, single-stranded RNA with 10 ORFs and several overlapping genes [8, 10, 11]. Despite its high structural and genetic similarity to SARS-CoV, the rapid spread and different clinical manifestations have raised new questions about the evolution of this new virus.

In South America (SA), the virus arrived at approximately the end of February 2020 from the European continent, followed by multiple additional introductions [12-16]. Despite control measures by different governments to reduce infections, the virus spread across the continent and, as was expected, new variants emerged due to mutation, increasing viral genetic diversity.

Genetic diversity is a relevant factor for viruses to adapt to changing environments [17]. In the case of infectious pathogens, DNA/RNA sequence analysis can be useful for understanding what factors (mutation, natural selection, genetic drift, recombination, migration [18, 19]) shape the genetic diversity patterns in natural populations, and this a key tool for understanding species evolution, monitoring the emergence of new variants, associating viral genetic variations with disease, and assessing the spread of an outbreak [19, 20]. Considering that population-genetics-based studies integrating sequence data from SA are limited, here, we analyze the diversity of SARS-CoV-2 from seven SA countries using methods based on genealogy, population genetics, and coalescent theory.

Materials and methods

Data collection

A total of 2098 SARS-CoV-2 genome sequences from SA (Argentina, 38; Brazil, 1070; Chile, 183; Colombia, 175; Ecuador, 90; Peru, 384; Surinam, 71; Uruguay, 75; Venezuela, 12) submitted by December 2020 were downloaded from the database of the Global Initiative on Sharing Avian Influenza Data (GISAID, https://www.gisaid.org/, Online Resource 1). These sequences were screened, and only those with less than 1% of the data missing were used. Additionally, virus genome reference sequences of the Wuhan SARS-CoV-2 (GenBank accession number NC_045512.2), bat BtCoV RaTG13 (China National Center for Bioinformation accession number GWHABKP00000001), bat BtCoV RmYN02 (GISAID accession number EPI_ISL_412977) and pangolin MP789 (GenBank accession number MT121216.1) were obtained.

Sequence alignments and analysis of genetic diversity

Only entries in which missing data or ambiguous nucleotides accounted for less than 1% of the genome sequence were included in the dataset. Ultimately, 1141 genome sequences (Argentina, 33; Brazil, 426; Chile, 162; Colombia, 107; Ecuador, 67; Peru, 303; Suriname, 43) were aligned with the Wuhan SARS-CoV-2, bat BtCoV RaTG13, BtCoV RmYN02, and pangolin MP789 coronavirus genome reference sequences using MUSCLE [21], followed by manual editing. The 5′ and 3′ ends of the aligned sequences were excised, resulting in an alignment consisting of 29,404 sites from ORF1a to ORF10. Genome sequences from Uruguay and Venezuela were omitted because they were submitted as partial genome sequences or had more than 1% of the data missing.

Several measures of DNA sequence polymorphism within or between populations were estimated, including the number of segregating sites (Ss), the number of singleton (S) and parsimony-informative sites (Ps), the number of haplotypes (H, the set of single-nucleotide polymorphisms found in the genome), the nucleotide polymorphism per site (θW), the nucleotide diversity per site (π) applying the Jukes-Cantor correction [22], and the average number of nucleotide differences between any two sequences (κ). A sliding window for nucleotide diversity with a window length of 100 sites and a step size of 25 sites was used. All of these measurements were made with DnaSP v.6 software [23]; gaps or sites with missing data were excluded. In addition, the Phylogenetic Assignment of Named Global Outbreak Lineages (PANGOLIN) tool v.3.1.5 [24] was used to identify the SARS-CoV-2 lineages in the data set.

Assessing molecular evolution

Natural selection, genetic drift, mutation, recombination, and gene flow are the essential explanatory factors used for understanding the population genetic basis of evolutionary change [18, 19]. Several tests were performed to assess these evolutionary forces. Tajima’s D [25] and Fu and Li’s D*/F* [26] frequency spectrum tests were performed to assess the null hypothesis that all mutations are selectively neutral [27]. In addition, haplotype-based tests, such as K, Hd [28], and Fu’s Fs [29], were also implemented. Statistical significance for the frequency spectrum and the haplotype-based test was evaluated using coalescence simulations, obtaining confidence intervals. All of these procedures were carried out using DnaSP v.6 software. Tajima’s D and Fu’s Fs were also computed on Arlequin v.3.1 software [30].

Natural selection at the inter- and intraspecies level was assessed using the coding alignment. First, the molecular evolution rate was inferred by partitioning the alignment at codon positions 1, 2, and 3 to obtain the relative substitution rates at all three codon positions using a Bayesian approach in BEAST v.1.7 software [31]. Then, the selection signatures (positive and purifying) at individual sites (codons) were detected using the single-likelihood ancestor counting (SLAC) [32], mixed-effects model of evolution (MEME) [33], and fast, unconstrained Bayesian approximation (FUBAR) methods [34]. A sliding window for omega rates (ω = dN-dS) was also used with the SLAC data. Because recombination can bias the inference of positive selection sites [35, 36], the genetic algorithm for recombination detection (GARD) method [37] was used before running the aforementioned tests. Later, the effective number of codons (ENC) [38] and the codon bias index (CBI) [39] were estimated to measure the selective pressure at the translational level [40, 41].

On the other hand, linkage disequilibrium (LD) was assessed by using the ZNS statistic [42]; this is the average of R2 over all pairwise comparisons. Linear regression between LD and nucleotide distance was then performed to ascertain whether recombination could have occurred. Recombination was also evaluated using the GARD algorithm [37] and RDP v.5 software [43].

Demographic history in South America

Statistically significant values in the neutrality test could be due to demographic factors. Therefore, the demographic history of SARS-CoV-2 in SA was assessed by computing the distribution of the observed number of differences between pairs of haplotypes (the mismatch distribution) using Arlequin software. The expected values for a demographic expansion were calculated and plotted against the observed values. Both the pure demographic expansion and the spatial expansion models were used to compute the sum of square deviations (SSD) between the observed and expected mismatch as a test statistic, as well as the raggedness index of the observed distribution defined by Harpending [44]. In addition, a Bayesian skyline plot (BSP) was used to examine changes in the effective population size (Ne) through time using the Bayesian Markov chain Monte Carlo (MCMC) approach available in BEAST v.1.7 software. BEAUti v.1.7 was used to generate a BEAST XML file for BSP using the sample collection dates, the generalized time-reversible (GTR) substitution model with gamma site heterogeneity (+G), base frequency estimates (+F), and a strict clock. The MCMC length chain was run for 50 million generations. Tracer v1.5 [45] was used check the effective sample size (ESS) as well as to construct the BSP plot for each population.

Genetic structure indices based on the allelic content of haplotypes and their frequencies [46] were estimated by molecular variance analysis (AMOVA). The pairwise FST index was computed for all population pairs; this can be used as short-term genetic distances between populations [47, 48]. We also tested the hypothesis of individuals’ random distribution between population pairs by performing an exact test of population differentiation [49]. These analyses were performed using Arlequin software. NETWORK v.10 software was then used for constructing a median-joining network [50] to evaluate possible mutational pathways giving rise to SARS-CoV-2 haplotypes (the set of single-nucleotide polymorphisms found on the genome), their distribution, and frequencies. The network was inferred with the star contraction [51] and post-processing MP options [52] to delete all superfluous median vectors and links not contained in the shortest trees in the network. In addition, the GTR+G+F model was used to infer a maximum-likelihood (ML) phylogenetic tree using the IQ-TREE web server [53]. Topology reliability was evaluated by bootstrap analysis, using 1,000 replicates.

Diversity of the spike, E, and M proteins

The viral spike, E, and M proteins are exposed to the immune system and are therefore the main proteins of interest for vaccine development. Thus, the diversity of these proteins at the amino acid level was assessed independently, since different substitutions within them may alter their antigenicity, which is a potential immune escape mechanism.

Results

SARS-CoV-2 genetic diversity in SA

A total of 1141 genome sequences from seven countries in SA were analyzed to evaluate the genetic diversity of SARS-CoV-2. In SA, there were 1260 segregating sites (Ss), 679 haplotypes (H), and 57 lineages (L). Of the Ss, 65.6% were singleton sites and 34.4% were parsimony-informative sites. The average number of nucleotide differences (κ = 8) was like previous estimates around the world [54, 55]. The diversity estimators displayed low values (θW = 5.7 × 10-3 and π = 2.8 × 10-4), with Peru, Argentina, and Suriname having the viral populations with the highest nucleotide sequence diversity; Colombia had the lowest diversity (Table 1). Because the sample collection dates differed from one country to another, these estimators were calculated again using only the sequences obtained between March and May 2020, and the results were similar to those reported above. The nucleotide diversity sliding window showed that the most diverse regions of the genome were located within the ORF1a, ORF3, E, ORF6, and N regions (Fig. 1).

Table 1

Estimators of SARS-CoV-2 genetic diversity in South America

	n	Sites	Ss	S	Ps	M	H	L	θw	π	κ	Tajima	Fu and L		Fu
												D	D*	F*	Fs
South America
	1142	28,945	1260	827	433	1278	679	57	0.0057	0.00028	8.10	− 2.78	− 19.2	− 10.22	< − 30.0
Estimators per country
Argentina	33	29,109	68	41	27	68	28	7	0.0006	0.00032	9.31	− 1.65	− 2.61*	− 2.51*	− 15.17
Brazil	426	29,369	547	420	127	550	300	14	0.0028	0.00022	6.46	− 2.82+	− 15.31+	− 9.85+	− 32.28
Chile	162	29,397	153	104	49	153	79	16	0.0009	0.00020	5.90	− 2.50+	− 7.82*	− 6.40*	− 80.93+
Colombia	107	29,314	121	85	36	121	58	15	0.0008	0.00015	4.40	− 2.66+	− 6.62*	− 5.90*	− 33.15+
Ecuador	67	29,409	71	40	31	71	26	11	0.0005	0.00024	7.06	− 1.78+	− 3.57+	− 3.44+	− 5.202
Peru	303	29,393	522	318	204	524	186	23	0.0028	0.00039	11.46	− 2.69+	− 9.87*	− 7.16*	− 24.63+
Suriname	43	29,409	70	39	31	70	35	2	0.0006	0.00030	8.82	− 1.60	− 2.71*	− 2.75*	− 21.97++

n number of sequences analyzed, Sites number of sites analyzed, Ss number of polymorphic segregating sites, S singleton sites, Ps parsimony-informative sites, M number of mutations, H number of haplotypes, L number of lineages, θw nucleotide polymorphism, π nucleotide diversity, κ average number of nucleotide differences

*p < 0.05, **p < 0.02, +p < 0.01, ++p < 0.001

Fig. 1

Sliding window of nucleotide diversity. SARS-CoV-2 diversity in seven South American countries was assessed by computing the nucleotide diversity per site. Window length, 100 bp; step size, 25 bp. A SARS-CoV-2 genome model is given below the sliding window displaying the ORFs.

Regarding lineages, the earliest A and B lineages were present in some countries. A was found in Chile, while B was present in Brazil, Chile, Colombia, and Peru. The lineages A.1, A.2, and A.5 were found in Chile, Colombia, and Peru, whereas the B.1 and B.1.1 lineages were found in all of the countries (Online Resource 2, Supplementary Table S1). The country with the most lineages was Peru, while Suriname was the country with the fewest lineages (Table 1 and Online Resource 2, Supplementary Table S1).

Several tests were performed to determine the impact of evolutionary forces on SARS-CoV-2 diversity in SA. Tajima’s D and Fu and Li’s D*/F* neutrality tests showed statistically significant negative values (Table 1, Online Resource 2, Supplementary Table S2). The haplotype number and haplotype diversity did not show significant differences from those expected under neutrality. On the other hand, Fu’s Fs estimator had statistically significant negative values. All of these estimators were then calculated locus by locus. Tajima’s D, Fu and Li’s D*/F*, and Fu’s Fs showed statistically significant negative values in almost all loci. Moreover, the haplotype number was higher than expected under neutrality at some of them (Online Resource 2, Supplementary Table S3).

Comparison of the relative substitution rates at the three codon positions showed that the third position had a much higher evolution rate (about 12-fold faster) than the other two codon positions (Online Resource 3). Furthermore, codon-based methods showed that more than 2660 codons were affected by purifying selection (Fig. 2 and Online Resource 4). However, positive selection signals were also found throughout the genome (Fig. 2 and Online Resource 4). Most of the positively selected codons were found in the regions encoding the non-structural proteins ORF1a (22 sites), ORF1b (9 sites), and ORF3 (9 sites), as well as in the region encoding the spike protein (9 sites). Most of these changes were non-conservative (i.e., the amino acid was replaced by one with different physicochemical properties), with D614G in the spike protein, G204R in N, and V/L3606F in ORF1 being the most frequent polymorphisms. Moreover, S84L was detected in ORF8, which is a reversion to the ancestral state. Finally, seven codons under positive selection showed substitutions corresponding to interspecies changes, all of which were non-conservative substitutions. These were in the regions that encode the spike (3 sites), ORF1a (1 site), M (1 site), ORF8 (1 sites), and N (1 sites) proteins.

Fig. 2

Sliding window for the ω (dN - dS) rate. The nonsynonymous (dN) and synonymous (dS) substitution rates throughout the SARS-CoV-2 genome were determined based on SLAC data. A SARS-CoV-2 genome model is given below the sliding window displaying the 10 ORFs. Values greater than zero indicate a positive selection, while values below zero indicate purifying selection.

The ENC and CBI were estimated to assess whether the low genetic diversity in SARS-CoV-2 could be (in part) due to selective pressure at the translational level. The average ENC value was 45.4, while the CBI value was 0.38 in all countries. Finally, some estimators were calculated to determine if recombination had taken place in the SARS-CoV-2 genomes in SA. The ZNS values were relatively low (Online Resource 2, Supplementary Table S2). Additionally, no correlation was observed between the LD and nucleotide distance (Online Resource 5A). Using the GARD method as well as the programs in RDP, no intraspecies recombination signals were found. However, RDP software did find recombination signals at the interspecies level (Online Resource 5B).

The demographic models based on the distribution of the observed number of differences between pairs of haplotypes in the sample displayed a unimodal distribution for each country (Figs. 3 and 4). The population data from Argentina, Brazil, and Peru showed that the observed distribution fit the expectations under the pure demographic expansion model (Fig. 3). Meanwhile, the data from Chile and Colombia fit the spatial expansion model better (Fig. 4). In Suriname, the observed data showed a multimodal distribution (Figs. 3 and 4).

Fig. 3

Mismatch distribution for the pure demographic expansion model. Frequency distributions of the observed number of pairwise nucleotide differences between haplotypes for seven South American countries are shown. The solid line is the theoretical distribution under the assumption of pure demographic expansion.

Fig. 4

Mismatch distribution for the spatial expansion model. Frequency distributions of the observed number of pairwise nucleotide differences between haplotypes for seven South American countries are shown. The solid line is the theoretical distribution under the assumption of spatial population expansion.

The effective virus population size (Ne) changes were then inferred within each country (Fig. 5). An increase in Ne was observed throughout January and February 2020; Brazil and Peru showed a drastic increase. The lowest Ne was observed in Suriname. In that country, Ne showed an increase around May 15. This behavior seems to agree with World Health Organization (WHO) statistics on the average number of cases reported daily for each country.

Fig. 5

Bayesian skyline plot (BSP) of the SARS‐CoV‐2 outbreak. BSP assessed changes in effective population size (Ne) in seven South American countries. The y‐axis indicates Ne, and the x‐axis shows the time in days/months. The solid black line represents the estimated median value, and the blue shades indicate the 95% highest posterior density.

In order to assess the population structure, an AMOVA analysis was performed evaluating the significance of the covariance components associated with the different possible levels of genetic structure. Given that the sequences recovered from the GISAID database were collected at different times, those collected during the same period for all countries (March to May) were chosen to perform this analysis. Both within and among viral populations, genetic structure indices showed statically significant values (Table 2). Eighty percent of the estimated variation in SARS-CoV-2 in South America was due to variation within populations, while differences between them could explain about 20% of the variation. The degree of genetic differentiation between populations was inferred using the FST index. Comparisons between populations showed statistically significant values higher than 0.05, indicating a slight genetic differentiation (Table 2). The exact test of population differentiation also displayed statistical support for genetic differentiation in SA. This also was evident when the lineage data were considered. There were several non-shared lineages circulating within each country (Online Resource 2, Supplementary Table S1).

Table 2

Analysis of molecular variance (AMOVA) analysis and population genetic differentiation estimated from March to May 2020 using the fixation index (FST)

	Source of variation	Sum of squares	Variance components	Percentage variation
	Among populations	503.17	0.73	19.73++
	Within populations	2717.86	2.97	80.27++
	Total	3221.03	3.70

Arg Argentina, Bra Brazil, Chi Chile, Col Colombia, Ecu Ecuador, Per Peru

**p < 0.02; +p < 0.01; ++p < 0.00001.

Likewise, the network results showed three high-frequency haplotypes (Online Resource 6A). The first one included the Wuhan-Hu-1 haplotype and sequences belonging to the B, B.1, B.1.8, B.1.9, B.1.110, B.1.111, B.1.67, B.1.195, B.1.223, B.1.225, B.1.547, B.3, and B.6 lineages. The second high-frequency haplotype included sequences from the B.1, B11, B.1.1.28, B.1.1.33, and N.4 lineages. The last high-frequency haplotype included sequences from the B.1, B.1.1, B.1.1.10, B.1.1.28, B.1.1.33, B.1.1.231, B.1.1.388, and B.59 lineages. Others derived from these three major haplotypes formed a star-like shape in the network. The low representation of haplotypes from Suriname is due to the fact that only 12 cases were reported from March to May, and only three sequences were available in the GISAID database. By including all of the sequences recovered by December 2020, a new haplotype network was inferred (Online Resource 6B) that preserved the topology shown in Online Resource 6A. It was also evidenced that all of the Suriname haplotypes (which cluster into the B.1.219 and N.2 lineages) seem to be exclusive to that country, showing a clear structuring. The ML phylogenetic tree agreed with these results (Online Resource 6C).

SARS-CoV-2 has four structural proteins that are essential for viral assembly, attachment, and entry into host cells. The spike protein is a 1273-amino-acid (aa) glycoprotein that is involved in the interaction with human receptors. One hundred fourteen haplotypes (“alleles”) of this protein were observed in South America, 109 of which had the D614G mutation (Online Resource 7). Haplotype 1 (Wuhan-Hu-1) had a low frequency, while haplotype 2 (with the 614G mutation) within the B.1 lineage was the most frequent in all populations except in Suriname, where haplotype 109 (E281V and D614G) within the B.1.219 lineage was the major haplotype.

Unlike the spike protein, M (222 aa) and E (75 aa) showed fewer haplotypes, with six and three, respectively (Online Resource 8). The Wuhan-Hu-1 haplotype was the most frequent in all countries for both the M and E proteins. However, haplotype 4 of M (D3G) reached 20.9% in Ecuador and 12.9% in Peru, and haplotype 2 of E (V58F) reached 18.2% in Argentina.

Discussion

Natural selection, genetic drift, and demographic processes play an important role in shaping the genetic background of natural populations, and the genetic diversity within them is a determining factor for survival of the species. Hence, characterization and monitoring of polymorphisms are key to understanding the dynamics and evolution of pathogens [19, 20]. The emergence of SARS-CoV-2 in humans is a relatively recent event that has triggered the current pandemic. The virus was reported in SA on February 26, when the first case of COVID-19 was registered in a 61-year-old traveler arriving in Brazil from Italy. A few days later, more cases of COVID-19 were reported in other countries in the region. Since then, it has spread to all SA countries.

A high diversity of lineages was found in SA, with B.1 (20.6%), B.1.1.33 (19.3%), B.1.1 (13.5%), and B.1.1.28 (11.4%) being the most frequent lineages. However, the genetic diversity could be higher than the lineage diversity, since some mutations might be present in some but not all sequences within lineage. For instance, in the P.2 lineage, the spike protein substitutions are E484K, D614G, V1176F, and F565L, but the last of these substitutions has been detected in some sequences but not all. Thus, SARS-CoV-2 genetic diversity was investigated. The first reports show low diversity and a low frequency of recurrent mutations [56-58]. According to our analysis, singleton sites (S, rare variants carried by a unique sequence) were 1.9-fold higher than parsimony-informative sites (Ps, recurring variants in two or more sequences), contrasting with the first reports, where S was three times higher than Ps [56-58]. This suggests an increase in the frequency of some variants that emerged after the appearance of this new betacoronavirus.

Various estimators can be calculated from the number of segregating sites (Ss) and the number of sequences to assess the genetic diversity of populations. Although nucleotide polymorphism (or Watterson estimator θW) was considered the best estimator of SARS-CoV-2 diversity [57], the nucleotide diversity (π) is a better estimator, since it considers the Ss frequency, and it is not influenced by the sample size. The nucleotide diversity in SA (π = 0.00028) was lower than in other continents/countries (Africa, π = 0.00047; Oceania, π = 0.00046; United States of America, π = 0.00044; Europe, π = 0.00043; Asia, π = 0.00042 [59]). All of these values were lower than those observed for other human viruses (H1N1 = 0.004-0.012 [60, 61]; H3N2 = 0.004 [60]; hepatitis B virus = 0.015 [62]; parainfluenza virus = 0.002 [63]), which could be due to the fact that SARS-CoV-2 is a new human virus.

The nucleotide diversity sliding window showed that the most diverse regions were located at the ORF1a, ORF3, and N loci (Fig. 1), agreeing with previous reports [59]. The viral populations from Peru, Argentina, and Suriname had the highest nucleotide diversity in SA. The covid performance (https://interactives.lowyinstitute.org/features/covid-performance/) gauges the relative performance of countries at different points in the pandemic. No correlation between nucleotide diversity (or lineage diversity) and covid performance was observed. Argentina (π = 0.0003), Brazil (π = 0.0002), and Colombia (π = 0.0001) got low scores (<10), while Peru (π = 0.0004), Ecuador, and Chile (π = 0.0002) had a score >30. The diversity observed could have resulted from different numbers of introductions [12-16] and the subsequent mutations emerging after arrival of the virus in each country, some of which have led to new lineages and variants of epidemiological interest [64].

The genetic diversity observed in natural populations results from five evolutionary forces: mutation, genetic drift, natural selection, recombination, and migration [18, 19]. Several tests were carried out to determine the impact of these forces on SARS-CoV-2 diversity in SA by the first months of the pandemic. Neutrally tests showed statistically significant negative values (Table 1). This could be the outcome of directional selection, selective sweep, or a population expansion. Selection is expected to affect diversity at some regions within loci. In contrast, a demographic process is expected to affect diversity throughout the genome.

This latter scenario is supported by the distribution of polymorphisms as well as by locus-by-locus analysis. The polymorphisms distribution in each population (the mismatch distribution) displayed a unimodal pattern (Figs. 3 and 4); populations that have experienced demographic growth in the recent past show this behavior [65]. Furthermore, statistically significant negative values for Fu’s Fs by country and by locus (Table 1 and Online Resource 2, Supplementary Table S3) indicate an excess in the number of recent mutations, supporting a demographic expansion. This is expected, since COVID-19 is a pandemic. However, Suriname showed a multimodal distribution, which could be due to either a population in equilibrium, co-circulation of multiple lineage, migration, or a population admixture.

Population expansion is expected to occur rapidly after the arrival of the virus. To address this scenario, changes in the effective population size (Ne) was inferred (Fig. 5). Ne increased before the first reported cases (February 26–March 13), suggesting that the virus may have arrived at SA in early January/mid-February. Phylogenetic analysis conducted in Colombia supports this [14, 16]. Once again, Suriname showed a different behavior. There, an increase in Ne seems to have occurred in mid-May although the first case was reported on March 13. Thus, Ne increase does not correspond to the arrival of the virus but to the control measures taken by the local government. They adopted social distancing and isolation only one day after the first case. This response appears to have contained virus dispersion. Both the Ne and case number increase coincide with the parliamentary elections, when citizens went out to vote on May 25 [66].

In contrast to Suriname, a drastic Ne increase was observed in Ecuador, Peru, or Brazil. From the time of the arrival of the virus in Brazil, the central government was discouraging the measures of physical distancing and lockdown brought in by state governors, and this, together with the complicated social situation, led to an increase in the number of cases [67], triggering an increase in the viral Ne. The aforesaid may suggest the effectiveness of early isolation to contain viral spread, and this should be considered in future outbreaks or new potential pandemics.

On the other hand, a population structuring is probable, considering the spatial expansion model observed in SA and since founder effects and genetic drift appear to be the main forces in SARS-CoV-2 evolution. Both within and among viral populations, genetic structure indices showed statically significant values (Table 2). The degree of genetic differentiation between populations using the FST index and population differentiation exact test showed low but statistically significant values, suggesting a slight genetic differentiation. Network and ML tree results and the presence of several non-shared lineages support this (Online Resource 2, Supplementary Table S1 and Online Resource 6).

Purifying selection also modulates SARS-CoV-2 evolution. The high evolutionary rate observed in the third codon position and the large number of sites with a dN - dS value < 0 confirm this (Fig. 2 and Online Resources 3 and 4). However, positive selection (+Sel) cannot be completely ruled out. Codon-based methods identified sites under +Sel at the intraspecies level. Most of them were found in regions encoding nonstructural proteins (Online Resource 4). Furthermore, the amino acid replacements change physicochemical properties. It is not clear why these proteins are the most diverse and had +Sel signatures. Thus, these variants should be monitored. Likewise, it needs to be evaluated through functional assays whether these replacements influence viral processes. Additionally, seven non-conservative amino acid replacements among SARS-CoV-2, bat CoV, and pangolin CoV displayed +Sel signatures (Online Resource 4). These could be ancient adaptive changes of a progenitor of SARS-CoV-2, as has been suggested previously [68].

Analysis of the structural spike protein displayed 114 unique haplotypes (“alleles”) in SA (Online Resource 7), and some sites were under +Sel, which could be due to diversifying selection, as the immune system exerts selective pressure or directional selection as an adaptation to the host. Within this last scenario, position D614G stands out, which was identified as being under positive selection in all of the populations analyzed here. Spike haplotypes containing both variants were found in SA, with the G mutation being the most frequent. The 614G polymorphism appears to affect viral infectivity and transmissibility [69-71]. Considering that a) under drift, polymorphism is expected to be fixed or get lost randomly, b) this variant had positive selection signals, and c) 614G has almost completely replaced the 614D variant worldwide, it is feasible to think that this mutation has been fixed by directional selection.

Recently, several SARS-CoV-2 variants and/or lineages have emerged and are circulating worldwide. Some of them are considered Variants of Interest (VOI: Epsilon [B.1.427 and B.1.429], Eta [B.1.525], Iota [B.1.526], Kappa [B.1.617.1], Zeta [P.2]), Variants of Concern (VOC: Alpha [B.1.1.7], Beta [B.1.351], Delta [B.1.617.2], Gamma [P.1]) or Variants of Note (A.23.1 and B.1.525). Fifty-seven different lineages from the dataset analyzed were found in SA, but none of these variants were detected. However, several spike mutations defining the Gamma and Zeta variants (i.e., E484K, F565L, D614G, and V1176F of Zeta and P26S, E484K, and D614G of Gamma) were already circulating in SA before the first report of these variants (Online Resource 7). Thus, virus isolates (or lineages) containing some of these mutations were their ancestors. For instance, the Gamma variant is a sublineage of B.1.1.28, which has been observed to be circulating in Brazil since March 2020. Additionally, other mutations (i.e., L5F, T95I, Q677H, P681R, and K1191N) observed in other VOIs or VOCs were present in the dataset analyzed here (Online Resource 7), suggesting a potential in SA for the emergence of new VOI or VOC, since they have been defined by multiple substitutions that usually arise independently within different lineages (convergent mutation) [72].

On the other hand, selection at the transcriptional level can also modulate diversity [40, 41]. To assess this scenario, ENC and CBI values were computed. ENC values ranges from 61 (all of the synonymous codons for each amino acid are used equally) to 20 (only one codon per amino acid). Likewise, CBI values range from 0 (uniform use of synonymous codons) to 1 (maximum codon bias). Since average values were higher than 45 and lower than 0.4, respectively, selective pressure at the transcriptional level was ruled out. Finally, recombination is an evolutionary force that is important for many viruses [73, 74]. Recombination signals have been observed previously in the SARS-CoV-2 genome [57, 75]. In addition, phylogenetic evidence of recombination has been reported when novel bat coronaviruses closely related to SARS-CoV-2 were considered [76]. Here, recombination events were observed at the interspecies level (Online Resource 5B), but not within the SARS-CoV-2 genome. No correlation was observed between the LD and nucleotide distance (Online Resource 5A). Moreover, GARD and RDP did not show intraspecies recombination signals. Thus, although recombination has taken place in coronaviruses, this mechanism is not yet a relevant factor in the evolution of the SARS-CoV-2 genome in SA.

SARS-CoV-2 has four structural proteins that are necessary for virus assembly, receptor binding, and entry into host cells, making them potential vaccine candidates. However, mutations within these proteins could alter their antigenicity. The sequences of the M and E proteins were found to be highly conserved in SA, and therefore, they might be considered good vaccine candidates. However, they appear to be poorly immunogenic [77], and their use as vaccine candidates is therefore limited. On the other hand, the spike glycoprotein, which plays an important role in the interaction with the ACE2 receptor, is the main target of the immune system and has therefore been used in vaccine design [77]. The vaccines developed so far use the full-length spike sequence or just the RBD domain. Mutations in RDB (D354, Y364, and F367) have been reported in several variants from Asia, Europe, and the USA, suggesting that these could impact antigenicity [62]. Mutational scanning analysis has detected several escape mutations within the N-terminal spike domain and RBD at positions 148, 150-151, 406, 409, 439, 443-452, 455-456, 472-473, 477, 483-486, 490, 493, 499, and 503 [72]. In SA, some of these sites (or sites close to them) as well as other sites throughout the protein showed polymorphism (Online Resource 7). It is not yet clear whether these variations in SA strains impact antigenicity. If so, the Sputnik and Janssen vaccines could be highly effective in SA, as RBD had a low degree of polymorphism among the sequences analyzed here. Conversely, the long-term effectiveness of the Pfizer, Moderna, and AstraZeneca vaccines, which induce antibodies directed against different regions of the spike [77], could be affected by the diversity found in the spike protein. Further studies are necessary to evaluate whether spike polymorphism will undergo a balancing selection process, maintaining an excess of variants at intermediate frequencies in the population as an immune evasion mechanism or due to selective pressure exerted by vaccination.

Below is the link to the electronic supplementary material.

Online Resource 1 Acknowledgement data source (XLSX 28 KB)

Online Resource 2 Supplementary Tables S1–S3 (PDF 165 KB)

Online Resource 3 The relative substitution rates at all three codon positions (TIF 911 KB)

Online Resource 4 Positions in the SARS-CoV-2 genome under natural selection in South America. The codons under positive or purifying (negative) selection were identified by codon-based methods for each population independently. The physicochemical properties of the amino acid replacements were assessed using the PAM250 matrix (Dayhoff and Schwartz, 1978). *, interspecies; +, reversion to ancestral character; italics, conservative change; bold, non-conservative change (XLSX 139 KB)

Online Resource 5 Recombination analysis. A. Relationship between linkage disequilibrium (LD) and nucleotide distance estimated for each South America country. The red line represents the regression line between LD and nucleotide distance. B. Recombination events found at the interspecies level (TIF 1563 KB)

Online Resource 6 Median-joining network for South America haplotypes and ML tree. A. Median-joining network for South America haplotypes. The figure shows the SARS-CoV-2 haplotypes identified in seven South American countries between March and May 2020. B. Median-joining network for South America haplotypes. The figure shows the SARS-CoV-2 haplotypes identified in seven South America countries using all sequences available by December 2020 in the GISAID database. For both networks, the star contraction option was chosen for simplifying interpretation of the network. Each node is a haplotype, and its size indicate its frequency. The lines connecting the haplotypes represent the different mutational paths explaining the relationship and evolutionary origin. The early A and B lineages are shown. C. SARS-CoV-2 phylogenetic tree inferred by maximum likelihood using genome sequences from South America. The topology was obtained using the GTR + G + F model. Lineage data are provided (TIF 2850 KB)

Online Resource 7 Spike haplotype alignment. The alignment shows the haplotypes (“alleles”) found in the spike protein. The numbers at the top indicate the amino acid positions where segregating sites were observed. The numbering is based on the SARS-CoV-2 Wuhan reference genome sequence (GenBank accession number NC_045512.2). The dots indicate identical amino acids. Red highlighted sites indicate the mutations observed in Variants of Interest or Variants of Concern (TIF 2931 KB)

Online Resource 8 M and E haplotype alignment. The alignment shows the haplotypes (“alleles”) found in M (A) and E (B) proteins. The numbers at the top indicate the amino acid positions where segregating sites were observed. The numbering is based on the SARS-CoV-2 Wuhan reference genome sequence (GenBank accession number NC_045512.2). The dots indicate identical amino acids (TIF 1133 KB)