Exploring the Interaction between E484K and N501Y Substitutions of SARS-CoV-2 in Shaping the Transmission Advantage of COVID-19 in Brazil: A Modeling Study.

The COVID-19 pandemic poses serious threats to global health, and the emerging mutation in SARS-CoV-2 genomes is one of the major challenges of disease control. Considering the growth of epidemic curve and the circulating SARS-CoV-2 variants in Brazil, the role of locally prevalent E484K and N501Y substitutions in contributing to the epidemiological outcomes is of public health interest for investigation. We developed a likelihood-based statistical framework to reconstruct reproduction numbers, estimate transmission advantage associated with different SARS-CoV-2 variants regarding the marking (identifying) 484K and 501Y substitutions (including Alpha, Zeta, and Gamma variants) in Brazil, and explored the interactive effects of genetic activities on transmission advantage marked by these two mutations. We found a significant transmission advantage associated with the 484K/501Y variants (including P.1 or Gamma variants), which increased the infectivity significantly by 23%. In contrast and by comparison to Gamma variants, E484K or N501Y (including Alpha or Zeta variants) substitution alone appeared less likely to secure a concrete transmission advantage in Brazil. Our finding indicates that the combined impact of genetic activities on transmission advantage marked by 484K/501Y outperforms their independent contributions in Brazil, which implies an interactive effect in shaping the increase in the infectivity of COVID-19. Future studies are needed to investigate the mechanisms of how E484K and N501Y mutations and the complex genetic mutation activities marked by them in SARS-CoV-2 affect the transmissibility of COVID-19.

INTRODUCTION

COVID-19, the etiological agent of which is severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), posed a serious threat to global health and swept the world in 2020; the pandemic is still ongoing., As of March 31, 2021, more than 127 million COVID-19 cases had been confirmed worldwide, with more than 2.7 million associated deaths.

The control of COVID-19 requires knowledge of the driving factors that may affect the transmission process,; virus mutation is one of the major challenges., Around September 2020, genetic variants carrying the N501Y substitution on the spike (S) protein of SARS-CoV-2 were first detected in the United Kingdom and, then spread globally and trended to reach fixation rapidly in many places (e.g., South Africa, Brazil, the United States, and the United Kingdom)., In Brazil, the 501Y variants, as well as other amino acid changes, were clustered into the B.1.1.28.1 lineage by COVID-19 Genomics Consortium UK, which is also known as the variant of concern (VoC) 202101/02. The B.1.1.28.1 lineage is a descendant of the B.1.1.28 lineage that has another similar descendant lineage, B.1.1.28.2, carries the E484K substitution but not N501Y substitution., The mutation E484K was first identified in South Africa and became prevalent in many places, including the United Kingdome and Brazil. These emerging variants may affect the epidemiological characteristics of COVID-19, and the protective effects of vaccines.– Considering the growing patterns of the epidemic curve in Brazil, the possible contributions of both E484K and N501Y substitutions are of public health interest for investigation.

The rapid spread of 501Y and 484K variants indicates a possible transmission advantage over their preceding variants. On one hand, recent analyses reported evidence that the N501Y substitution was associated with an increase in infectivity of COVID-19,,,– which appears similar to the situation of D614G substitution reported previously.– On the other hand, the relationship between E484K substitution and COVID-19 transmissibility appears inconclusive., The survival or functional profile of pathogen could be altered through genetic mutation and, as a consequence, change its infectivity. Referring to the previous studies on seasonal influenza viruses, a few key amino acid substitutions may lead to changes in antigenic features and epidemiological outcomes,, and the interaction among them may become more complicated. As an example, the R384G substitution in the nucleoprotein (NP) of H3N2 virus enhances the ability of in-host immune escape, which increases transmissibility, but this substitution appears detrimental. In contrast, the commutations including E375G and M239V in NP could compensate and restore the viral fitness or functionality of H3N2 virus,, just as the mutated strains rapidly reached fixation in 1993–1994 influenza season. For the COVID-19 epidemics in Brazil, how the mutation activities marked by E484K and N501Y substitutions, as well as the possible interactive effect between them, might shape the transmission advantage remains unassessed.

Exploring the role of mutation activities in determining disease transmissibility is of importance to understand how the evolutionary process at molecular scale may shape the epidemiological outcomes at population scale.,,, In this study, we adopt a statistical framework to infer the real-time transmissibility associated with different SARS-CoV-2 variants with respect to E484K and N501Y substitutions in Brazil. We explore the interactive effects between E484K and N501Y in shaping the transmission advantage of COVID-19.

METHODS

SARS-CoV-2 sequencing data and COVID-19 surveillance data.

The SARS-CoV-2 strains were obtained via the global initiative on sharing all influenza data (GISAID) with collection dates ranging from January 1, 2020 to January 31, 2021, in the Brazil. A total of 4,210 complete human SARS-CoV-2 strains were retrieved. We excluded sequences with more than 5% ambiguous amino acids during the alignment, and a total of 4,052 sequences were included for further analysis. Multiple sequence alignment was performed using MAFFT version 7, and the “Wuhan-Hu-1” (GISAID: EPI_ISL_402125 or GenBank: NC_045512.2) SARS-CoV-2 genome is considered as the reference sequence.

The surveillance data of COVID-19 cases in the Brazil were collected via the WHO COVID-19 surveillance platform. To avoid the under-ascertainment due to reporting delay, we drop the observations since February 2021. As such, the surveillance data of COVID-19 cases from January 1, 2020 to January 31, 2021, are included in the analysis, which match the period of SARS-CoV-2 sequencing data. To adjust for the weekly cycle in the COVID-19 case time series, the 7-day moving average is adopted for further analysis. The COVID-19 cases time series are shown in Figure 1A.

Figure 1.

(A) The daily number of COVID-19 cases in Brazil, (B) the reconstructed reproduction number (R), (C) proportions of the 484E/501N variants, (D) 484E/501Y variants (including B.1.1.7 or Alpha variants), (E) 484K/501N variants (including P.2 or Zeta variants), and (F) 484K/501Y variants (including P.1 or Gamma variants). Panel B shows the estimated Rs of 484K/501Y variants (in blue) and the other three types of variants (in gray). The dots are the estimates, and bars are the 95% confidence intervals (CIs). In panels C through F, the dots are the observations, the curve indicates the mean fitting results, and the shading area indicates the 95% CI. This figure appears in color at www.ajtmh.org.

Statistical parameterization.

Variant-specific reproduction number.

The time-varying reproduction number is commonly adopted to quantify the instantaneous transmissibility of infectious disease in an epidemic. Using the estimation framework in Cori et al., the epidemic growth is modeled as a branching process. Thus, the reproduction number at time t, R(t), is expressed as the ratio of C(t) over , which is commonly known as the renewable equation., Here, the C(t) is the observed new incidences of COVID-19 at time t. The function w(·) is the distribution of the generation time (GT) of the disease, that is, COVID-19. The GT is defined as the time interval between the time of exposure, that is, being infected, of a primary case and that of his associated secondary case in the consecutive transmission generation. Thus, the distribution w(·) is predefined in our model, which is commonly estimated from contract tracing surveillance data.– To set up the analysis for COVID-19, we consider w as a Gamma distribution having mean (±SD) values of 5.3 (±2.1) days by averaging the GT estimates from the existing literatures.– Slight variation in the settings of the GT will not affect our main findings.

To incorporate the information of SARS-CoV-2 variants, we denote the proportion (or prevalence) of the j-th variant of concern (VoC) at calendar time t by ρ(t), which is time-varying. Straightforwardly, for all t. We denote the variant-specific reproduction number for the j-th VoC at time t by R(t), and we have . We consider the expected new incidences infected by j-th VoC at time t, denoted by E[C(t)], which can be modeled in Eq. (1).

Equation (1) will be used to formulate the likelihood function in the remaining sections. Note that the index j is merely used for labeling instead of ranking.

Transmission advantages and their interactive effects.

For convenience, we label the original variant as the (j =) 0-th VoC, that is, the 484E/501N variant, and thus its associated variant-specific reproduction number is R = 0. Similarly, we label the 484E/501Y, 484K/501N, and 484K/501Y variants as the j = 1, 2, and 3 VoC, respectively.

Following previous work,, the transmission advantage of the mutated variant against the original type is defined as the ratio (denoted by η) of the strain-specified reproduction numbers. Considering the transmissibility of the original 484E/501N variant (R = 0) as the reference level, the transmission advantage of the j-th VoC is . Thus, the reproduction number of cases infected by the j-th VoC is R = η·R=0, and η = 0 = 1 by definition. If η > 1, the j-th variant may be more infectious than the original genetic variant, and vice versa. In addition, the overall reproduction number is .

We consider η as a constant, which reflects an intrinsic nature of the j-th SARS-CoV-2 variant, and thus η is invariant with time. Hence, we have R(t) = η·R=0(t) for all time t. Then, we calculate the expected number of COVID-19 cases, E[C(t)], at time t in Eq. (2).

As such, for the observed sequencing data, the expected chance (or probability) that a randomly selected strain at the t-th day is j-th VoC, E[ρ(t)], is given in Eq. (3).

Straightforwardly, for all t.

For the E484K and N501Y substitutions, the 484E/501Y (j = 1, i.e., including the B.1.1.7 or Alpha variants) and 484K/501N (j = 2, i.e., including the P.2 or Zeta variants) are two variants with merely one substitution, whereas the 484K/501Y variant (j = 3, i.e., including the P.1 or Gamma variants) has both substitutions. To explore the interactive effects on the variant-specific reproduction number, we compare η=3 and the product of (η=1·η=2). If η=3 is larger than (η=1·η=2), the E484K and N501Y substitutions may enhance the transmissibility than their separated partial effects, and vice versa.

Likelihood-based inference.

According to Eq. (2), we construct the likelihood function of the daily number of cases using a Poisson-distributed framework with observation at C and rate parameter at E[C] as in Eq. (4).

Here, the C is the observed number of COVID-19 cases on day t and is the discretized C(t), which means . The value of C can be obtained from the number of COVID-19 cases time series as shown in Figure 1A. Note that the superscript (c) merely indicates that the likelihood function is for the number of cases, which does not indicate the power.

For the observed sequencing data, we denote the numbers of j = 0, 1, 2, and 3 variants by m=0,, m=1,, m=2,, and m=3,, respectively, for the day t. Thus, we model the sampling process of genetic variants using a generalized Bernoulli distribution, that is, categorical distribution, with probabilities at E[ρ(t)]s in Eq. (3). The likelihood function is constructed in Eq. (5).

Here, the E[ρ,] is the expectation of variant prevalence for day t. Note that the superscript (s) merely indicates that the likelihood function is for genetic variants, which does not indicate the power.

With Eqs. (4) and (5), we reconstruct the R=0, time series, denoted by , and estimate η=1, η=2, and η=3 using the overall log-likelihood function ℓ defined in Eq. (6).

We calculate the maximum likelihood estimates (MLE) of parameters to determine transmission advantage of 484E/501Y (η=1), 484K/501N (η=2), and 484K/501Y (η=3) variants by using the likelihood framework defined in Eq. (6). The 95% confidence intervals (95% CI) are calculated using the profile likelihood estimation framework with a χ2 quantile as the cutoff,, which has also been adopted in previous work.–

Fitting schemes and their selection.

We explore the interactive effects between E484K and N501Y substitutions of SARS-CoV-2 in shaping the transmission advantage in Brazil. We consider the following eight fitting scenarios with respect to η and investigate the role of E484K and N501Y substitutions contributing to the transmissibility of COVID-19:

• all of η=1, η=2, and η=3 are assumed at 1:

scenario (#1): η=1 = η=2 = η=3 = 1;

• two of η=1, η=2, or η=3 are assumed at 1, and the remaining one is freely estimated:

scenario (#2): η=1 = η=2 = 1,

scenario (#3): η=1 = η=3 = 1, and

scenario (#4): η=2 = η=3 = 1;

• one of η=1, η=2, or η=3 is assumed at 1, and the remaining two are freely estimated:

scenario (#5): η=1 = 1,

scenario (#6): η=2 = 1, and

scenario (#7): η=3 = 1;

• none of η=1, η=2, or η=3 is assumed at 1, and all three of them are freely estimated:

scenario (#8): none of η=1, η=2, or η=3 is assumed at 1.

The settings of the eight fitting scenarios are presented in Table 1.

Table 1

The summary of transmission advantage estimates of different types of SARS-CoV-2 variants in Brazil under different scenarios

Scenario	Transmission advantage of				AIC	AICc	BIC	HQIC	Remarks
	Original variant	New emerging variant
	484E/501N	484E/501Y	484K/501N	484K/501Y
(#1)	1 (reference)	1 (assumed)	1 (assumed)	1 (assumed)	6895.6	6960.3	8915.2	7615.1	Baseline model
(#2)	1 (reference)	1 (assumed)	1 (assumed)	1.23 (1.04–1.41)	6888.7	6953.8	8914.4	7610.4	Two types of new variant are assumed to have no effect
(#3)	1 (reference)	1 (assumed)	1.05 (0.93–1.17)	1 (assumed)	6893.5	6958.6	8919.3	7615.2
(#4)	1 (reference)	1.08 (0.78–1.41)	1 (assumed)	1 (assumed)	6893.5	6958.6	8919.3	7615.2
(#5)	1 (reference)	1 (assumed)	1.12 (1.01–1.22)	1.28 (1.08–1.47)	6888.3	6953.8	8920.3	7612.2	One type of new variant is assumed to have no effect
(#6)	1 (reference)	1.15 (0.83–1.52)	1 (assumed)	1.21 (1.05–1.42)	6889.9	6955.4	8921.9	7613.8
(#7)	1 (reference)	1.08 (0.79–1.43)	1.04 (0.92–1.16)	1 (assumed)	6895.3	6960.8	8927.2	7619.2
(#8)	1 (reference)	1.26 (0.92–1.66)	1.14 (1.03–1.27)	1.33 (1.13–1.56)	6888.4	6954.3	8926.5	7614.5	Full model

AIC = Akaike information criterion; AICc = corrected AIC for small sample size; BIC = Bayesian information criterion; HQIC = Hannan-Quinn information criterion. The highlighted scenario (#2) is selected as the main result, and shown in Figure 1.

We conduct the model fitting and parameter estimation under each scenario. The scenario with the best fitting performance is selected according to lowest values of 4 different information criteria including Akaike information criterion (AIC), corrected AIC for small sample size (AICc), Bayesian information criterion (BIC), and Hannan-Quinn information criterion (HQIC).

Sensitivity analysis.

Sensitivity analysis was conducted to examine the robustness and significance of the determine transmission advantage estimates, that is, ηs. We examine the consistency of both directions of the effects and their 95% confidence intervals (CI) under alternative settings. The following three sensitivity checking schemes are performed.

For the first scheme, we consider a univariate logistic regression model between the overall R as response and ρ=1,, ρ=2,, and ρ=3, as regressors. The regression coefficients of all ρ, are evaluated as the effect size of SARS-CoV2 variants on the overall transmissibility of COVID-19. For the second scheme, we repeat the estimating process of transmission advantage with alternative PDF of GT, that is, w(·), which is introduced in the previous section, “Variant-specific reproduction number.” We consider shorter and longer versions of mean GT at 4 days,, and 7.5 days, respectively. For the third scheme, we repeat the analysis by replacing the Poisson-distributed likelihood function in Eq. (4) with a negative binomial (NB) distribution to further account for the superspreading potential of COVID-19 transmission. For the setting of NB distribution, we fixed the dispersion parameter at 0.4, which follows the estimation in recent studies.,–

RESULTS AND DISCUSSION

As of this writing, in Brazil, the epidemic curve had grown since March 2020 with two major epidemic waves, the first in August 2020, and the second is ongoing (Figure 1A). The original 484E/501N variants started being replaced by the other three types of variants in October 2020 and had almost vanished in Brazil after January 2021 (Figure 1C). In January 2021, the prevalence of the emerging 484K variants was 71.1% and 46.6% for 501Y variants. The prevalence of 484E/501Y (j = 1 type) variants is 10.1%, 34.6% for 484K/501N (j = 2 type) variants, and 36.5% for 484K/501Y (j = 3 type) variants (Figure 1D–F). As such, the linkage disequilibrium (LD) is calculated at 0.03, which indicates the occurrence of two mutations is likely random. Specially, the 484K/501Y variants are classified into the B.1.1.28.1 (or P.1, or 20J/501Y.V3) lineage.

For the eight fitting scenarios summarized in Table 1, we find that the transmission advantage of 484K/501Y (η=3) variants is estimated larger than 1 significantly and consistently. In contrast, the scale of η=1 for 484E/501Y variants appears statistically unclear compared with 1, that is, not significantly larger than 1. For the model selection, scenario 2 has the lowest AICc, BIC, and HQIC, and scenario 5 has the lowest AIC, and AICc. Both scenarios 2 and 5 also have close values of AIC with a difference of only 0.4. As such, scenario 2 is considered the main result and thus is presented in Figure 1B.

The modeling framework in this study links the mutation activity at molecular scale and COVID-19 transmissibility at population scale. We reconstruct the instantaneous reproduction numbers (R) of COVID-19 cases infected by 484K/501Y variants and other variants in Brazil under fitting scenario 2 (Figure 1B). The overall trends of reproduction numbers are relatively high in the early phase of outbreak before and in May 2020 and the November 2020 for the second major epidemic wave, but gradually decrease thereafter. The average scale of reproduction number during the early outbreak is largely consistent with previous estimates.,,–

We report the estimated proportions of four types of SARS-CoV-2 variants, E[ρ], fit the observed sequencing data well (Figure 1C–F). We infer the transmission advantage η=3 for 484K/501Y variants at 1.23 (95% CI: 1.04–1.41), which means the E484K and N501Y substitutions together increase 23% of COVID-19’s transmissibility in Brazil, whereas other emerging variants, 484E/501Y and 484K/501N, are unlikely to have significant transmission advantage. Thus, in Figure 1B, the reproduction number of the 484K/501Y variant appears higher than that of the other genotypes (non-484K/501Y, or non-Gamma variants). For sensitivity checking, we find that the η=3 estimates are consistently and significantly larger than 1 in similar scales as the main estimates (data not shown), which validates our findings.

We focus on the second major epidemic wave because E484K and N501Y substitutions emerged during the same period. Although the reproduction number of non-484K/501Y variants has fluctuated around 1 since December 2020, the reproduction number, R=3, of 484K/501Y variants are largely greater than 1 during the same period, which has led to a large epidemic wave in Brazil in 2021 (see Figure 1A). Given that 484K/501Y variants trend to reach fixation, both the herd immunity threshold and intrinsic growth rate of epidemic may increase, which was previously discussed regarding the situation in United Kingdom,, and thus the local nonpharmaceutical interventions for COVID-19 control may be enforced. Hence, we highlight the importance of our analytical framework, such that the public health risks related to viral mutations may be controllable with early preparedness.

The increase in transmissibility associated with the 484K/501Y variants is biologically reasonable. The N501Y substitution is a mutation on a key contact residue in the receptor binding domain on the S protein and is found to increase the ability of human angiotensin-converting enzyme 2 binding and cell infectivity in animal models, which appears similar to the previous D614G substitution. The E484K substitution is also located in the viral RBD and confers resistance to several monoclonal antibodies by affecting the binding process.,, However, according to the η=1 estimates in Table 1, N501Y substitution alone appears less likely to form an advantage in fitness without accompanying 484K. Similarly, E484K substitution alone also might not secure a concrete transmission advantage robustly or significantly. Our finding indicates that the combined impact of 484K/501Y outperforms their independent effects, which implies statistically an interactive relationship.

Previous studies reported a higher case fatality risk among the individuals infected by SARS-CoV-2 strains in the B.1.1.7 lineage,, some of which carry both E484K and N501Y substitutions. Given the transmission advantage of 484K/501Y (i.e., η=3 > 1), the increasing intensity of COVID-19–related mortality is a public health concern. Clinical severity remains largely unassessed for the B.1.1.28 lineage in Brazil, and unexpected clinical outcomes may warrant adjustments in the treatment strategies. The emergence of 484K/501Y variants and its mutations (e.g., K417N or V1176F), combined with other VoC in Brazil and other places, implies the capacity of SARS-CoV-2 to evolve new phenotypes rapidly. Although the neutralizing level of BNT162b2 vaccine-elicited sera is recently found to be satisfactory for 484K/501Y variants, further investigation is required for other vaccine candidates at population scale.

This study has the following limitations. First, our analysis was based on the sequence data released in GISAID and thus is subject to the selection bias of sequences being released to the public domine. Second, the reconstruction of reproduction numbers relies on the setting of the generation time (GT). Theoretically, the GT distribution might be altered by the mutated strains. However, by screening the literature, we find no evidence that GT is associated with the E484K or N501Y substitution in SARS-CoV-2, and thus we model GT distribution of COVID-19 w(·) as a fixed Gamma distribution, following previous studies.– Third, we consider w(·) as a fixed distribution. In the real-world situation, the time interval between transmission generations might vary,, which may affect the reconstruction of the reproduction number. However, the long-term trends of R estimates are unlikely to change due to slight variation in GT., Thus, we consider the impact of this limitation on the inference of transmission advantage may be negligible, and our model can be extended to a more complex context with the time-varying GT data available. Fourth, due to the lack of data from different Brazilian regions, we aggregated the national COVID-19 cases in Brazil to reconstruct the reproduction number series (Figure 1B). We acknowledge this analytical scheme neglects the heterogeneities in epidemiological characteristics of COVID-19 transmission,,,, geographic separation in SARS-CoV-2 variants, individual response and vulnerability to COVID-19,– and various nonpharmaceutical interventions, for wide regions across different Brazilian locations. We note that the transmission advantage may vary under different local settings or situations. Fifth, ideally, C(t) in the R estimation should be the number of COVID-19 cases with onset at time t. However, because surveillance data by date of onset are unavailable, we adapted the current dataset by reporting data as a proxy for the COVID-19 incidence time series. If one considers a constant reporting lag, the R estimates will have the same trends but are shifted for this lag. Considering that a similar reporting delay also occurred during the collection of SARS-CoV-2 sequencing data, the effects of the two reporting lags could be counteracted. We believe that this approximation is unlikely to affect the main conclusions of this study. Furthermore, with detailed reporting lag of information for each individual case, adjustments for reporting delay can be carried out based on our current analytical framework. Sixth, this study focuses on exploring the effects on changing disease transmissibility associated with mutation activities, but real-world biological mechanisms, which are usually more complex, remain uncovered. Future studies are needed to explore the mechanisms of how E484K and N501Y mutations in SARS-CoV-2 affect the transmissibility of COVID-19. Seventh, there exist other mutations in the B.1.1.28 lineage, such as K417N and V1176F, but we merely considered the E484K and N501Y because they are dominant in B.1.1.28.1 and B.1.1.28.2 lineages, respectively. Given the lack of individual patient information, time-series data were used in this work, which means there is information loss from the data aggregation. As pointed out in previous studies, the independent effects of each commutation may not be disentangled in this study due to identification issues that the samples might fail to inform each estimate. Eighth, the interpretation of our findings should be limited to the COVID-19 epidemics in Brazil, but similar investigations could be conducted for other regions. Ninth, for simplification, we consider the transmission advantage (η) of new variants versus the wildtype as constant over time. This model assumption may not necessarily be strictly held considering several real-world determinants, including the accumulation of population immunity against different strains, selection pressure due to intervention strategies,, and behavioral factors related to disease spread and transmission. Alternatively, the estimating framework of η can be extended into a real-time basis. Lastly, as a data-driven study, the estimated association should be treated with caution. In an ecological setting, although our analysis provides statistical evidence on the likelihood of causality, the findings in this study cannot guarantee causality, which requires further biomedical experiments for verification.

Supplemental materials