Full vaccination against COVID-19 suppresses SARS-CoV-2 delta variant and spike gene mutation frequencies and generates purifying selection pressure

Is there correlation between mutation frequency and vaccination?

To explore this question, we first analyzed the correlation between the rates of full vaccination and the point mutation frequency (Mf) of COVID-19 delta variants’ genome in 20 countries. Complete SARS-CoV-2 genome sequences with high coverage from June 20 to July 3 2021 in 20 countries: Australia (N=121), France (N=788), Germany (N=955), Indonesia (N=97), India (N=171), Ireland (N=617), Israel (N=333), Italy (N=642), Japan (N=105), Mexico (N=368), Netherland (N=456), Norway (N=142), Portugal (N=782), Singapore (N=131), Spain (N=689), Switzerland (N=131), Sweden (N=786), Turkey (N=428), United States (N=537), and United Kingdom (UK, N=3534) were collected from the Global Initiative on Sharing All Influenza Data (GISAID) (https://www.gisaid.org/). Mf was calculated as Pi / (Ln × Ns). Pi is the number of instances of polymorphism detected within the genome/locus. Ln is the nucleotide length of the genome/locus. Ns is the number of sequenced entities present in the dataset (Roy et al., 2020). The information of fully vaccinated rate on 3 July 2021 in individual counties were obtained in Our World in Data (https://ourworldindata.org/covid-vaccinations).

We found that Mf was logarithmically reduced as the full vaccination rate increased in 16 of 20 countries (R²=0.878, P<10⁻⁷, Figure 1A). To our knowledge, this is the first evidence suggesting that vaccinations could successfully suppress viral mutations. Since the spike protein is the target of the vaccination program, we further examined mutations of delta variant spike gene. Likewise, we found that the Mf of the spike gene was also logarithmically reduced as the full vaccination rate increased this time in 17 of the 20 countries (R²=0.829, P<10⁻⁶, Figure 1B). With a 10.8% vaccinated rate, Mf was exceptionally low in Australia, likely as a result of the strict state lockdown restrictions. In contrast, the Mfs are higher in Japan, Switzerland, and United States, suggesting that their social mitigation strategies have been less successful. This also suggests that he mutation frequency may be not only affected by vaccination only, but may also reflect overall control of the virus.

• Download figure
• Open in new tab

Figure 1.

Correlation between full vaccinated rate (https://ourworldindata.org/covid-vaccinations) and mutation frequency (Mf) of the full genome (A) or the spike gene (B) of SARS-CoV-2 delta variants from June 20 to July 3 2021 in 20 countries: Australia (AUS), France (FRA), Germany (GER), Indonesia (IDA), India (IND), Ireland (IRL), Israel (ISR), Italy (ITA), Japan (JPN), Mexico (MEX), Netherland (NED), Norway (NOR), Portugal (POR), Singapore (SGP), Spain (ESP), Switzerland (SUI), Sweden (SWE), Turkey (TUR), United States (USA), and UK. Logarithmic regression (solid) line was draw based on 16 countries (pink) with a calculated 95% confidence interval (dashed lines). Japan, Switzerland, USA, and Australia are labeled in different colors as outliers.

How does vaccination affect selection pressure?

SARS-CoV-2 is subject to many powerful selection pressures which may contribute to vaccine escape. Like other viruses, the delta variant appears to have been selected for its increased transmissibility regardless of vaccination. Previously we reported on the mutations and selection pressure of SARS-CoV-2 among passengers in close-quarters quarantine on the Diamond Princess cruise ship between January and March 2020 (Yeh and Contreras, 2021). Until then, it had been unclear whether COVID-19 variants were evolving randomly or under selection pressures generated by vaccination and/or other mitigation efforts.

To understand how selection pressure operates within COVID-19 virus populations, we analyzed the polymorphism of the full SARS-CoV-2 genome using neutrality tests based on the site-frequency spectrum: (1) Tajima’s D test, (2) normalized Wu and Fay’s DH test, and (3) Zeng’s E test, using DNASP6 software with the bat coronavirus RaTG13 as the outgroup sequence (Rozas, et al., 2017). The Tajima D test compares pairwise nucleotide diversity (π, the average number of nucleotide differences per site between two sequences) and total polymorphism to infer selection and/or demographic events (Tajima, 1989). DH test is sensitive to the changes in high-frequency variants and primarily affected by directional selection and no other driving forces, such as demographic expansion (Zeng et al., 2006). Zeng’s E test contrasts high and low frequency regions of frequency spectrum, so this test can detect population growth after a sweep (Zeng et al., 2006). Statistical significance of observed values of all tests was obtained by coalescent simulation with free recombination after 10,000 replicates. The probability for each test statistic was calculated as the frequency of replicates with a value lower than the observed statistic (two-tailed test) by DNASP6 (Rozas, et al., 2017).

We observed that Tajima’s D values were negative and significantly deviated from zero (−2.199 to -2.878, P < 0.05) among 19 countries, indicating an excess of variants of low frequency (Table 1). Because the negative Tajima’s D values could result from demographic expansion and/or purifying/positive selection, we further employed normalized DH (insensitive to population growth) and Zeng’s E test (Zeng et al., 2006). Because fixation causes quick loss of high-frequency alleles by drift, and excess of low-frequency mutations are generated, Zeng’s E test is very sensitive to population growth. Genomic polymorphism of these viral populations revealed significantly negative DH values (−1.177 to -2.791, P<0.05), but E values were not significantly different from zero (−0.04 to 0.482, P>0.05) (Table 1). Similar results were obtained after confining our analysis to the spike gene of SARS-CoV-2 delta variants using Tajima’ D, DH, and Zeng’s E tests (Table 1). Taken together, these data suggest that directional selection was operating in these countries, and demographic expansion was less likely the major driving force.

Selection pressure has been observed in unvaccinated SARS-CoV-2 populations, such as in the aforementioned Diamond Princess shipboard lockdown (Yeh and Contreras, 2021). Our current results lead to the hypothesis that vaccination could contribute to a part of the selection pressures on SARS-CoV-2, especially on the spike protein. It has been shown that, within populations, the frequency distribution of non-synonymous polymorphisms is negatively skewed relative to the distribution of synonymous polymorphisms under purifying selection (Hahn et al, 2002; Hughes, et al. 2005). Therefore, this hypothesis can be examined using a modified Tajima’s D statistics, which takes more negative values for non-synonymous (D_nonsyn) than for synonymous sites (D_syn) of spike gene under purifying selection (Hahn et al, 2002, Hughes, et al., 2005). The average number of pairwise synonymous differences (kS) and the average number of nonsynonymous pairwise differences (kN), the number of synonymous segregating sites (S_S), and the number of nonsynonymous segregating sites (S_N) were computed with equation as Austin Hughes described (1) S*_S=S_S/a1, (2) S*_N=S_N/a1, where a1 is as (Tajima, 1989, Hughes, et al., 2005). D_syn was defined as kS − S*_S, divided by the standard error of that difference, and D_nonsyn was defined as kN − S*_N, divided by the standard error of that difference (Tajima, 1989, Hughes, et al., 2005). The values of K_a/K_s (the ratio of the number of nonsynonymous substitutions per non-synonymous site (K_a), to the number of synonymous substitutions per synonymous site (K_s), in a given period of time), and the neutrality index of the McDonald-Kreitman (MK) test of the spike protein (RaTG13 as the outgroup sequence) were analyzed by DNASP6 software (Rozas, et al., 2017).

Unlike Tajima’s D test, it has been demonstrated that D_nonsyn and D_syn analysis is independent of sample size, so the value of D_nonsyn or D_syn can be compared between data sets of different sizes (Hughes et al., 2008). The same excess of low-frequency alleles in non-synonymous polymorphism was also observed (D_nonsyn = -1.482 to −2.714; P < 0.01). In addition, the K_a/K_s values were less than 1 and the neutrality index values of the MK test were significantly more than 1 (P<0.05) in all countries (Table 1). These results confirmed that purifying selection led to constraint on the available neutral mutations at non-synonymous sites of the spike gene of delta variants, consistent to previous reports of SARS-CoV-2 variants in 2020 (Chaw et al, 2020).

Because non-synonymous and synonymous mutations are affected unequally by purifying selection, D_nonsyn is expected to be disproportionally lower than D_syn (Hahn et al., 2002). It has been proposed that the value of ΔD (D_syn-D_nonsyn) increased as purifying selection became stronger, and ΔD (D_syn-D_nonsyn) is also able to rule out the homogenous effects (population expansion, selective sweep etc.) in neutrality tests (Hahn et al, 2002). We observed that ΔD of the spike gene of SARS-CoV-2 delta variants is positively proportional to fully vaccination coverage rate (R²= 0.505, P=0.001). When using D ratio (D_nonsyn/D_syn), the standard error terms in both equations of D_nonsyn and D_syn cancel out the sample size effect (Hughes, 2008). The positive correlation of D ratio and fully vaccination rate became even more evident (R²= 0.723, P<10⁻⁴) (Figure 2). Therefore, purifying selection pressure of SARS-CoV-2 spike gene increased as the vaccination coverage rate increased. This suggests that vaccination, in combination with other mitigation strategies, plays an important role in the purifying selection force of SARS-CoV-2 spike protein.

• Download figure
• Open in new tab

Figure 2.

Correlation between full vaccinated rate and D ratio (D_nonsyn /D_syn) of the spike gene of SARS-CoV-2 delta variants. Exponential regression (solid) line was draw based on 19 countries (pink square) with a calculated 95% confidence interval (dashed lines). Australia (AUS) is labeled as outliers.

Application of neutrality tests and their limitations

Neutrality tests are very powerful tools to study population genetics and viral evolution. However, these tests require careful interpretation due to two known limitations in this study: (1) sampling across countries and the representation of sequences in public databases, (2) the timing of introductions and vaccination campaigns. Previously we have applied various neutrality tests to analyze SARS-CoV-2 transmission and evolution in Diamond Princess cruise (Yeh and Contreras, 2021). In contrast to the shipboard quarantine, where the virus spread in a closed environment, many more complicated factors (mitigation efforts, travel, herd immunity, etc.) could influence viral evolution during pandemics among different countries.

Bhatt et al. have shown that the Tajima D test performs very well with low type I error rate (less than 5%) when the θ value (the substitutions per site of a sequence alignment) is between 0.1 to 10⁻⁴ in different RNA viruses (Bhatt et al., 2010). Because θ values in this study ranged between 0.0009 to 0.032, the type I error rate is low. The Tajima D test is also much less influenced by RNA recombination, which we have reported previously (Yeh and Contreras, 2020, Yeh and Contreras, 2021). Tajima D test results depend on the interplay of many parameters: positive and negative selection, population size dynamics, spatial structure and migration, random genetic drift, etc. A key limitation of the Tajima D test is the difficulty to separate the effects of each parameter. However, the selection pressure and demographic expansion can be verified in combination with other neutrality tests. Our results of DH and Zeng’s E test suggests that purifying selection, not demographic expansion, was the acting force. Tajima D may show negative values due to a recent bottleneck event rather than selection, but the bottleneck effect should impact all types of polymorphism equally (Tajima, 1989). Significant differences between D_nonsyn and D_syn (either ΔD or D ratio (D_nonsyn/D_syn)), which efficiently eliminated the homogenous effects, strongly indicates that purifying selection was at play and is positively corelated with vaccination coverage rate within these samples.

The inter-species divergence K_a/K_s test has been shown to be poorly equipped to detect positive selection in short time scales for the comparison of young evolutionary lineages like SARS-CoV-2 (Mugal, et al., 2014). Consistent with other analytical attempts, our K_a/K_s data was inconclusive (Table 1). This has been the case with other applications of the K_a/K_s test in SARS-CoV-2 spike gene, which resulted in extreme values and controversial results (Kang et al., 2021), In contrast to K_a/K_s, the MK test cannot be used to locate specific sites under selection, although this test enables us to estimate the adaptive substitution rate of a gene (Smith and Eyre-Walker, 2002). Variants of the MK test have been used to analyze the evolution of RNA viruses (Bhatt et al., 2010). Since the MK test relies on the comparison of the ratio of polymorphism to fixed differences of both synonymous and nonsynonymous mutations, selection of an outgroup sequence is critical for this test. The number of fixed synonymous mutations may be underestimated in a rapidly evolving virus like SARS-CoV-2 when a distant outgroup is applied (Baudry and Depaulis, 2003). Since the K_a/K_s or the MK tests may not be the appropriate tests for SARS-CoV-2 delta variants, we only included them to confirm the observation of purifying selection of SARS-CoV-2 spike gene here as well. We also did not include neutrality tests based on haplotype distribution or linkage disequilibrium, because they are expected to be strongly affected by RNA recombination that commonly occurred in SARS-CoV-2 populations (Yeh and Contreras, 2020, Ignatieva et al, 2021, Yeh and Contreras, 2021).

Perspectives

More virulent strains of SARS-CoV-2 have emerged with enhanced transmissibility and immune evasion properties. For example, multiple variants have escaped neutralizing antibodies developed to target the spike protein receptor-binding or N-terminal domain (Harvey, et al., 2021). What’s more, the case numbers of breakthrough infection caused by the delta variant have increased drastically worldwide (Farinholt et al., 2021).

Of all countries studied, our observations show that higher vaccination rates corresponded with lower mutation frequencies and higher D_nonsyn/D_syn values for the whole genome and the spike gene. We conclude then, that as the vaccination coverage rate increases, purifying selection forces of nonsynonymous mutations also increase.

Thus, we recommend that: 1) universal vaccination should be administered as soon as possible to suppress the generation of deadly mutations; 2) mitigations strategies such as personal protection equipment, social distancing, etc., should be continued until full vaccination is reached to prevent viral transmission; 3) genomic surveillance should be undertaken to monitor for new mutations; and 4) more sequence summary statistics of RNA virus evolution are required to facilitate understanding new COVID-19 outbreaks.