Estimating the strength of selection for new SARS-CoV-2 variants

Abstract

Controlling the SARS-CoV-2 pandemic becomes increasingly challenging as the virus adapts to human hosts through the continual emergence of more transmissible variants. Simply observing that a variant is increasing in frequency is relatively straightforward, but more sophisticated methodology is needed to determine whether a new variant is a global threat and the magnitude of its selective advantage. We present three methods for quantifying the strength of selection for new and emerging variants of SARS-CoV-2 relative to the background of contemporaneous variants. These methods range from a detailed model of dynamics within one country to a broad analysis across all countries, and they include alternative explanations such as migration and drift. We find evidence for strong selection favoring the D614G spike mutation and B.1.1.7 (Alpha), weaker selection favoring B.1.351 (Beta), and no advantage of R.1 after it spreads beyond Japan. Cutting back data to earlier time horizons reveals large uncertainty very soon after emergence, but that estimates of selection stabilize after several weeks. Our results also show substantial heterogeneity among countries, demonstrating the need for a truly global perspective on the molecular epidemiology of SARS-CoV-2.

Figure 1:

Profile likelihoods and credible intervals for the selection parameter s of the stochastic epidemic model for United Kingdom (UK), Netherlands (NL), and Japan (JP) and variants D614G, B.1.1.7, B.1.351, and R.1. The black dots indicate log-likelihood estimates of the fitted model with the corresponding fixed value of s, and the black curve is a smoothing spline through these log-likelihoods (see Methods). The dashed line shows the maximum-likelihood estimate, the gray box shows the 95% CI. The red intervals show the results of the population genetics model for these countries (cf. Fig. 5).

Figure 2:

Mechanistic model fit to D614G data in the Netherlands and UK. Panels A and B show many realizations of the prevalence of the background, I_wt, and D614G variant, I_mt for the maximum likelihood model fit. The insets show a close-up of the first weeks of the epidemics. Panels C and D show the number of deaths accumulated up to the week scale (black dots) and the model fits to those data. The bars around the points indicate the 95% predictive interval for the data according to the stochastic model. Panels E and F show the proportion of sequenced genomes with a glycine on position 614 of the spike protein in a given week. The blue lines are realizations of the model fit to the data. Vertical bars on the data indicate the 95% confidence intervals (CI) for the proportion based on the number of sampled genomes.

Figure 3:

Mechanistic model fit to B.1.1.7 data in the Netherlands and UK. Panels A-H are as in Fig. 2, but now for the B.1.1.7 variant. Panels G and H show the profile-likelihood results for the time-horizon analysis at a sequence of dates. The likelihood profiles (cf. Fig. 1) are plotted vertically, where one unit in log-likelihood space corresponds to a day in calendar time. The likelihood profiles intersect with the gray vertical lines at the boundaries of the 95% CI. The horizontal red line indicates s = 0. The CIs marked with a star (*) have a lower bound above s = 0.

Figure 4:

Mechanistic model fit to B.1.351 data in the Netherlands and R.1 in Japan. Panels A-F are as in Fig. 2 and panels G and H as in Fig. 3. Notice that the y-axes in panels E and F do not range from 0 to 1.

Figure 5:

Selection coefficients for each country from the population genetic model. Results are for the final time-points shown in Figs. 2–4. Points mark the median, and thick and thin lines are 50% and 95% CIs, respectively. Corresponding estimates of migration are in Fig. S2.

Figure 6:

Comparison of the population genetics model and the isotonic regression method. Each variant is analyzed at the final time horizon shown in Fig. 8, and results are shown for each country with sufficient data by that time.

Figure 7:

Estimated global distribution of selection coefficients for each of the four variants from the population genetics model. Our hierarchical model estimates the mean (top row) and standard deviation (bottom row) of a normal distribution from which the selection coefficient, s, of each country is drawn. In each panel, dark vertical lines mark the median, and the 90% CIs are shaded.

Figure 8:

Estimates of the global selection coefficient, from the hierarchical population genetics model, for each variant for differing amounts of data. (When only one country was present in the data, the non-hierarchical equivalent model was fit instead.) Dates show the first day on which that corresponding number of cases of that variant was first reached globally; they are the time horizon to which the data were cut back for that estimate. The number of days of data is measured from the time horizon back to a first day for each variant: 2020–01-04 for D614G, 2020-09-20 for B.1.1.7, 2020-10-01 for B.1.351, and 2020-10-24 for R.1. Points mark the median, and thick and thin lines are 50% and 95% CIs, respectively.