V367F Mutation in SARS-CoV-2 Spike RBD Emerging during the Early Transmission Phase Enhances Viral Infectivity through Increased Human ACE2 Receptor Binding Affinity

In total, 32,123 SARS-CoV-2 isolates with whole-genome sequences available in public databases and sampled before 31 March 2020 were analyzed. Among them, 302 isolates with RBD mutations were identified compared with reference strain Wuhan-Hu-1 that was the first isolated in Wuhan in December 2019. All the mutants were clustered into 96 mutant types, six of which were dominant mutant types that were found in more than ten isolates (Table 1): V483A (35×), V367F (34×), V341I (23×), N439K (16×), A344S (15×), and G476S (12×). V483A and V367F accounted for 11.59% and 11.26% of 302 mutants, respectively. V367F mutants occurred on 22 January 2020, which was the earliest dominant mutant type. The detailed alignments of amino acid sequences in the RBD of the mutants are shown in Fig. S1 in the supplemental material. These mutants were emergent in multiple continents, including Asia, Europe, North America, and Oceania. Most RBD mutants were circulating in Europe and North America (Fig. 1). V367F mutants were found on all four continents.

As SARS-CoV-2 S protein mediates attachment of the virus to cell surface receptors and fusion between virus and cell membranes, the polymorphism and divergence of the S gene were analyzed by DnaSP6 (version 6.12.03) (37). Overall, the S2 subunit was more diverse than the S1 subunit (Fig. 2A). The D614G mutation in the S1 subunit accounted for most of the polymorphism and divergence. Most of the other synonymous and nonsynonymous mutations occurred in the S2 subunit (Fig. 2B), of which the number of synonymous mutations was more than that of nonsynonymous mutations. During the early transmission phase, RBD was conserved and not as diverse as the S2 subunit (Fig. 2A and B). The six major mutant types are shown in Fig. 2C.

Since RBD is the only domain to bind the human ACE2 receptor and, in turn, initiates cell entry, it is believed that the RBD should be highly conserved. To test this hypothesis, we investigated the selective pressures on the S gene during the early transmission phrase by calculating the nonsynonymous/synonymous substitution rate ratios (dN/dS ratios) for various segments of the S gene from the 32,123 SARS-CoV-2 isolates (18, 19). The entire S gene exhibited a dN/dS of 0.86, close to the 0.83 reported previously (20) (Table 2). However, different regions of the S gene may be subject to different selection pressures. The S1 gene exhibited a high dN/dS (2.05) due to the widely spreading nonsynonymous mutation D614G in the S1 gene (21). The S2 gene exhibited a lower dN/dS (0.2388551), indicating that the S2 gene was more conserved than S1 under negative selection. The S1-RBD showed a lower dN/dS (0.7258567) than the entire S gene (0.8569338). Therefore, the functional relevance of the RBD mutations may be inferred.

To estimate the functional changes suggested by the RBD mutations, we performed molecular dynamics (MD) simulations for the prototype SARS-CoV-2 (Wuhan-Hu-1 strain) and the RBD mutants in order to assess their binding energy to the human ACE2 receptor. Each simulation was performed at 10 ns, and each model was simulated in triplicates. All trajectories reached a plateau of root mean square deviation (RMSD) after 2 to 5 ns (Fig. 3A), indicating that their structures reached equilibrium. All of the subsequent computations on thermodynamics were based on the 5- to 10-ns trajectories. The ΔG of the V367F mutant was significantly low (∼60 kJ/mol) (P = 0.0151), approximately 25% lower than for the prototype strain (−46.5 kJ/mol, calculated from the experimentally measured equilibrium dissociation constant [K_D]) (Fig. 3B), suggesting a significantly increased affinity to human ACE2; the other mutants showed a similar ΔG to that of the prototype (Fig. 3B). Compared to the K_D (14.7 nM) of the prototype RBD, the K_D of the V367F mutant was calculated as 0.11 nM (Fig. 3C), which was 2 orders of magnitude lower than for the prototype strain, indicating an increased affinity to the human ACE2 receptor. In comparison, the bat CoV RBD (strain RaTG13, with the highest genome similarity) showed a much lower binding affinity (K_D = 1.17 mM; ΔG = −17.4 kJ/mol) to human ACE2 than the pangolin CoV (K_D = 1.89 μM; ΔG = −33.9 kJ/mol). The affinity of the pangolin CoV to human ACE2 was lower than that of the SARS-CoV-2 prototype strain (K_D = 14.7 nM; ΔG = −46.5 kJ/mol) (Fig. 3B and C).

To elucidate the structural basis of the increased affinity of the V367F mutant, we investigated the dynamics of the residues in these structures in greater detail. The binding surface of the RBD to ACE2 is largely arrayed in a random coil conformation, which lacks structural rigidity. Logically, a firm scaffold should be necessary to maintain this conformation of the interaction surface and thus may facilitate the binding affinity. The beta-sheet structure scaffold, centered by residues 510 to 524 (Fig. 4A, in red), apparently provides this rigidity. “Higher-affinity” mutant V367F showed a considerable decrease of the root mean square of fluctuation (RMSF) at this region, demonstrating a more rigid structure; this was not observed for other mutants (Fig. 4B). Coincidentally, the substitutions that account for the affinity increase of V367F are all located near this fragment. Indeed, residues 475 to 485, which is a random coil near the binding site, showed a remarkably higher RMSF for the “similar-affinity” mutants than the higher-affinity V367F mutant (Fig. 4B). Moreover, the higher-affinity mutant V367F exhibited a generally decreased ΔG in the binding site region in contrast to that of the similar-affinity mutants (Fig. 4C).

As of 31 March 2020, among all of the mutants with dominant mutations in spike RBD (>10 isolates), most mutations resulted in a substitution of amino acids with similar properties. V367F is the only mutant with a higher binding affinity, as calculated by MD simulation (Fig. 3). Therefore, the binding affinity and the infectivity of the V367F mutant were further validated experimentally. First, we performed experiments to assess the binding affinity in vitro with a receptor-ligand binding enzyme-linked immunosorbent assay (ELISA) using purified S proteins and human ACE2 protein. The result showed that the V367F mutation lowered the 50% effective dose (ED₅₀) concentration (ED₅₀ = 0.8 ± 0.04 μg/ml) compared to that of the prototype (ED₅₀ = 1.7 ± 0.14 μg/ml) (Fig. 5A). This demonstrates that the V367F mutant has a higher affinity to human ACE2 than the prototype. Second, we performed surface plasmon resonance (SPR) experiments, which yielded the same conclusion: the prototype had a K_D of 5.08 nM compared to the V367F mutant with a K_D of 2.70 nM (Fig. 5B). Additionally, we performed a virus-cell interaction experiment to investigate the invasion efficiency of S proteins using an HIV backbone pseudovirus assay. ACE2-overexpressing Vero and Caco-2 cells were infected by the pseudoviruses bearing the mutant or wild-type RBD with the same multiplicity of infection (MOI). A higher infection efficiency is represented by the increased copy number of the integrated lentivirus genome. At 24 h postinfection (h p.i.), the V367F pseudoviruses showed 6.08× higher copy numbers than the prototype in Caco-2 cells (P < 0.01). At 48 h p.i., the V367F pseudoviruses showed 6.61× and 9.16× higher copy numbers than the prototype in Vero (P < 0.0001) and Caco-2 cells (P < 0.0001), respectively (Fig. 5C). Therefore, the computation and protein and cell validations were consistent with each other: the V367F mutant had enhanced affinity and infectivity.

The D614G mutation is located in the S1 region and is outside the RBD of the SARS-CoV-2 spike protein. It has been confirmed that the D614G mutants, which have spread widely, increase virus infectivity by elevating its sensitivity to protease (21). The V367F mutants were initially discovered in January 2020 in Hong Kong. Afterwards, the V367F mutants emerged mainly in Europe, including the United Kingdom, the Netherlands, Austria, and Iceland as well as in the United States, Australia, and China. The D614G+V367F dual mutant initially emerged in March 2020 in the Netherlands (see Table S1).

The phylogenetic analysis of the V367F mutant genomes during the early transmission phase showed that V367F mutants clustered more closely with the SARS-CoV-2 prototype strain in the L clade. Intriguingly, all of the dual-mutation variants (V367F+D614G) emerging later formed a distinct subcluster in GH, GR, and G clades, separate from the L and S clades (22) (Fig. 6). Compared with all the single-mutation (V367F) variants, dual-mutation (V367F+D614G) variants were emergent in both GH and GR clades. Multiple dual-mutation (V367F+D614G) variants located in different clades had mutations undetected in early V367F mutants. This indicates that the dual mutations may have evolved through several individual evolution events (Fig. 6).

Full-length gene sequences of SARS-CoV-2 were downloaded from the NCBI GenBank database, and the GISAID EpiFlu database (http://www.GISAID.org). In total, 34,702 SARS-CoV-2 full-genome sequences isolated during early transmission phase (samples collected before 31 March 2020) were downloaded, and the sequences with amino acid mutations in the S protein and RBD region were parsed. The genome sequences with either the V367F mutation or the V367F/D614G dual mutations in the S protein RBD were screened and analyzed in this study (see Table S1 in the supplemental material). For evolution analysis, 1,753 full-genome sequences from each NextStrain (26) clade were filtered and cluster random sampled from the GISAID EpiFlu database.

Alignment of S protein sequences from different sources and comparison of ACE2 proteins among different species were performed using MAFFT version 7 with default parameters (https://mafft.cbrc.jp/alignment/server/) and BioEdit (27, 28). Selection pressure analyses were conducted using the Nei-Gojobori method (Jukes-Cantor model) with Mega X (version 10.0.2) (18, 19). Substitution mutation analyses of mutants were performed and compared with the Wuhan-Hu-01 sequence using BioAider (version 1.314) (29). Phylogenetic analyses were conducted with IQ-Tree2, applying the maximum likelihood method with 1,000 bootstrap replicates using the GTR+F+R3 model by MFP ModelFinder (30). The trees were annotated and modified using Figtree (version 1.4.4) (https://github.com/rambaut/figtree/releases).

The complex structure of the SARS-CoV-2 S-protein RBD domain and human ACE2 was obtained from the Nation Microbiology Data Center (identifier [ID] NMDCS0000001) (PDB ID 6LZG) (31). Mutated amino acids of the SARS-CoV-2 RBD mutants were directly replaced in the model, and the bat/pangolin CoV RBD domain was modeled using SWISS-MODEL (32). Molecular dynamics simulation was performed using GROMACS 2019 with the following options and parameters: explicit solvent model, system temperature of 37°C, all-atoms optimized potentials for liquid simulations (OPLS/AA) force field, and LINCS restraints. With 2-fs steps, each simulation was performed at 10 ns, and each model was simulated three times to generate three independent trajectory replications. Binding free energy (ΔG) was calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method (GitHub; https://github.com/Jerkwin/gmxtool), with the trajectories after structural equilibrium assessed using root mean square deviation (RMSD) (35, 36). To calculate the equilibrium dissociation constant (K_D) and ΔG, the formula ΔG = RTlnK_D was used. Estimated ΔGs of the RBD mutants were normalized using the ΔG of the prototype strain, which was derived from experimental data.

The SARS-CoV-2 prototype S gene lacking the furin cleavage site, TM, and cytoplasm domains was cloned into pNPM5 vector (Novoprotein, NJ, USA) and fused with a C-terminal His₆ tag. The mutation was introduced by site-directed mutagenesis using the ClonExpress MultiS One Step cloning kit (Vazyme) and confirmed by PCR and sequencing. The wild-type and V367F constructs were transfected into HEK293 cells using polyethyleneimine. Since the S protein included a signal peptide in its N-terminal 14 amino acids, it was secreted into the medium. The expressed S proteins were purified from filtered cell supernatants with a Ni-nitrilotriacetic acid (NTA) column. Eluted protein solution was then dialyzed against phosphate-buffered saline (PBS; pH 7.4) for subsequent assays.

Human ACE2 protein was immobilized onto a microtiter plate at 5 μg/ml (100 μl/well). Each S protein (prototype and V367F) was added as a ligand at different concentrations, ranging from 0.03 μg/ml to 10 μg/ml, and then incubated for 2 h at 37°C to allow receptor-ligand interaction. The ligand-receptor mixture was then washed three times. One hundred microliters of horseradish peroxidase (HRP) anti-His tag antibody (BioLegend, USA) (diluted 1:20,000) was added to each well and allowed to react for 1 h. After three washes, the signal was visualized using 3,3′,5,5′-tetramethylbenzidine (TMB) solution (Sigma-Aldrich, USA), and a microtiter plate reader recorded the optical density at 450 nm (OD₄₅₀).

The SPR experiments were performed in a Biacore T200 instrument (GE, USA). For this, the SARS-CoV-2 S-proteins, either prototype or V367F, were immobilized on the sensor chip NTA (GE, USA), according to the manufacturer’s protocol. Human ACE2 protein was injected in each experiment at seven concentrations (6.25, 12.5, 25, 50, 100, 200, and 400 nM). For each cycle, the absorption phase lasted for 120 s and the dissociation phase lasted for 600 s. After each cycle, the chip was regenerated by washing with 350 mM EDTA and 50 mM NaOH for 120 s so that the chip was ready for the next round of S protein immobilization. Blank controls with 0 nM ACE2 were performed, with the blank signals subtracted from the cycle signals. All experiments were performed at 37°C. K_D values were calculated by fitting the curves using the software provided with the instrument.

The full-length S gene of SARS-CoV-2 strain Wuhan-Hu-1 (NC_045512.2) was cloned into a SARS-CoV-2 spike vector (PackGene, Guangzhou, China) and confirmed by DNA sequencing. Plasmid SARS-CoV-2 spike (G1099T), incorporating the V367F mutation in the S gene, was constructed by site-directed mutagenesis using the ClonExpress MultiS One Step cloning kit (Vazyme), as per the manufacturer’s protocol.

Generation of SARS-CoV-2 S HIV backbone pseudovirus was performed as previously described with some modifications (33, 34). Briefly, 293T cells, at approximately 70% to 90% confluence, were cotransfected with 9 μg of the transfer vector (pLv-CMV), 6 μg packaging plasmid (psPAX-lentiviral), and 6 μg envelope plasmid (pCB-spike or pCB-spikeV367F). Pseudoviruses were harvested at 48 h posttransfection and subsequently filtered into either Falcon or microcentrifuge tubes via a syringe-driven 0.45-μm filter. The harvested pseudoviruses were treated with recombinant DNase I (RNase-free) (TaKaRa, Japan) at 37°C for 10 min to eliminate carryover of plasmid DNA. Afterwards, the DNase I in the virus stock was inactivated at 95°C for 10 min after DNA digestion. Virus titration was measured by reverse transcription-quantitative PCR (RT-qPCR) targeting the WPRE gene of pseudoviruses using the Hifair III One Step RT-qPCR SYBR green kit (Yeasen), as per the manufacturer’s protocol. After reverse transcription (10 min at 42°C) and initial denaturation (5 min at 95°C), PCR amplification was performed for 40 cycles (15 s at 95°C, 60 s at 60°C). Primers were as follows: WPRE-F, 5′-CACCACCTGTCAGCTCCTTT-3′; WPRE-R, 5′-ACGGAATTGTCAGTGCCCAA-3′.

To detect S variant-mediated viral entry, Vero E6 and Caco-2 cells (5 × 10⁴) grown in 24-well plates were infected with either S-V367 or S-F367-bearing purified pseudovirus (pretreated with DNase I) (1 × 10⁶ pseudovirus RNA copies). At 24 h and 48 h postinfection, the cells were washed three times, and the total DNA was extracted. Relative fold increase of infected virus titers was calculated according to the WPRE DNA copy number of the lentiviral proviruses measured by TB Green Premix Ex Taq II real-time PCR kit (TaKaRa). Data from all of the samples were obtained from three independent experiments, and each sample was tested in triplicates.