shared this story
SARS-CoV-2 is a once-in-a-century pandemic, having emerged suddenly as a highly infectious viral pathogen. Previous phylogenetic analyses show its closest known evolutionary relative to be a virus isolated from bats (RaTG13), with a common assumption that SARS-CoV-2 evolved from a zoonotic ancestor via recent genetic changes (likely in the Spike protein receptor binding domain – or RBD) that enabled it to infect humans. We used detailed phylogenetic analysis, ancestral sequence reconstruction, and molecular dynamics simulations to examine the Spike-RBD’s functional evolution, finding to our surprise that it has likely possessed high affinity for human cell targets since at least 2013.
Viral pathogens are a continuous and evolving challenge for human populations.1,2 Most known viruses maintain species-specific infectivity, often co-evolving with their host to mirror animal species trees.3,4 While less common, the emergence of novel viral pathogens is of particular interest because they often exhibit abnormal degrees of infectivity and/or virulence,5 having not evolved to a natural selection balance with their new host.6 Viruses of animal origin include periodic Ebola outbreaks,7 the 1918 “Spanish Flu”,8 and most recently, SARS-CoV-2, the viral agent that causes COVID-19.9 In these cases, viruses spread through human populations after evolving to “cross the species barrier”.10 Yet, many questions remain for viruses of non-human origin: How do they acquire the ability to infect humans? Is it wholly dependent on “recognition” (a function typically mediated by protein-protein binding between viral capsid and target host cell), or must there be changes in other viral replication mechanisms as well? And specifically focusing on SARS-Cov-2, did it evolve to infect humans via many evolutionary changes or only a few? Was it dependent on a key functional shift in its ability to bind human cells, or is there evidence that other genomic changes were needed for it to acquire its strikingly high degree of infectivity? Answering these questions is critical if we are to understand both the origin of specific viruses, such as SARS-CoV-2, as well as the capacity of animal viruses to evolve human infectivity.
SARS-CoV-2 emerged in late 201911 and has high infectivity, spreading rapidly around the world, causing a global health emergency.12 A member of the Coronaviridae family of polymorphic, enveloped, single stranded RNA viruses,13 it is thought that SARS-CoV-2 evolved from a zoonotic origin,14,15 owing to its clear evolutionary relationship with coronaviruses that have been isolated from animals16 (its closest known evolutionary relative is the bat coronavirus, RaTG1317–20 and the second-closest known relative is a pangolin coronavirus, Pangolin-CoV).21 While most of the SARS-CoV-2 genome is most similar to the RaTG13 genome, some genomic regions, including the Spike glycoprotein Receptor Binding Domain (RBD) (which mediates viral entry into host cells), have greater sequence similarity to the Pangolin-CoV homolog,22 prompting some to suggest SARS-CoV-2 may be the product of recombination during co-infection.21–24
The Spike protein is a key component of the SARS-CoV-2 infection pathway.25 Knockout and overexpression studies have demonstrated that binding of the Spike-RBD to human angiotensin converting enzyme 2 (hACE2) mediates cellular entry of SARS-CoV-2.26–30 The protein sequence of this surface receptor is variable, with particular rare variants increasing patient susceptibility to SARS-CoV-2 infection.31 The SARS-CoV-2 Spike protein has been shown to bind the hACE2 receptor with greater affinity than the SARS-CoV-1 homolog, which has been suggested as a possible explanation for its greater infectivity.29 Additionally, many other related coronaviruses have been shown to be unable to bind hACE2 with sufficient affinity to support infection, raising the possibility that high hACE2 is a recently acquired trait for SARS-CoV-2.32–34 Given this, a critical question remains to be answered: How and when did the SARS-CoV-2 Spike protein evolve its relatively higher affinity for the hACE2?
With this question in mind, we set out to rob
ustly characterize the evolutionary changes that accompanied the emergence of SARS-CoV-2, and that distinguish it from its closest zoonotic relatives. We focused on the evolution of the Spike-RBD by leveraging its known evolutionary relationships to other animal and human viruses and employed ancestral sequence reconstruction in conjunction with molecular dynamics simulations to identify biochemical and biophysical changes that enhanced Spike binding to the hACE2 receptor.
Our initial phylogenetic analysis utilized whole viral genomic data, and generally supports prior studies’ conclusions, finding similar levels of nucleotide identity to the RaTG13 genome (96.0% sequence identity) and the Pangolin-CoV genome (90.0% sequence identity) (Supplementary Figure 1A).21,29 We next quantified the degree of evolutionary diversification that has occurred during SARS-CoV-2’s global spread. We performed an in-depth analysis of 479 sequences collected between December 30, 2019 and March 20, 2020, and observe 16 polymorphisms, including 11 missense mutations present in >5% of infections (Supplemental Table 1), each mapping to unique phylogenetic branches (Figure 1A). One monophyletic clade was primarily isolated within the United States and Canada, and is defined by two synapomorphic missense mutations: c.17848A>G and c.28134C>T.35 Since these occur in one of the most variable parts of the coronavirus genome, it is likely that its distribution is due to a founder effect and that it does not confer an evolutionary advantage.36 It is also worth noting that neither appears in the Spike protein-coding region, making it unlikely to impact hACE2 affinity.
We subsequently sought to investigate the evolution of SARS-CoV-2 infectivity by performing ancestral sequence reconstruction for the Spike-RBD amino acid sequence (Figure 1B). While cross-species protein sequence comparisons have been used to investigate critical amino acid changes in the SARS-CoV-2 Spike protein,37 leveraging the phylogenetic relationships allows us to infer the most likely ancestral sequences for this protein allows us to focus on the subset of genetic changes that are specific to its recent evolution.38 We inferred statistically well-supported reconstructions of the Spike-RBD sequence for both the common ancestor of all human SARS-CoV-2 cases (labelled “N1”, Figure 1B,D) and the its common ancestor with the closest animal virus (labelled “N0”, Figure 1B,D). N1 is identical to the Spike-RBD in the SARS-CoV-2 reference sequence, as expected, while the N0 Spike-RBD sequence is, to our knowledge, unique, reflecting the uniqueness of SARS-CoV-2’s viral origin.21,39 N0 differs from N1 at 4 positions (346, 372, 498, and 519 – Figure 1C). The reconstruction of N1 for each of those positions is statistically well-supported, with a posterior probability (P.P.) of 1 obtained from two independent calculations (Supplemental Table 2; Supplementary Methods). The reconstruction for N0 has high statistical support for positions 346, 372, and 519 (P.P. > 0.94), while position 498 was ambiguously reconstructed, with two alternate states comparably probable (Supplemental Table 2). All other positions were reconstructed with high confidence (P.P.>0.9). Together, these four changes (t346R, t372A, h/y498Q, and n519H) differentiate the evolved SARS-Cov-2 Spike protein from the most recent common ancestor with animal viruses (Figure 1B). As such, this ancestral virus must have existed at least as early as 2013 (as one of its descendants – RaTG13 – was isolated in that year), meaning that the branch between the N0 and N1 ancestors covers at least 7 years of molecular evolution (Figure 1B).
To quantify functional differences between the N0 ancestor and the Spike-RBD sequences, we conducted 10 ns molecular dynamics simulations (see Supplementary Methods) of the Spike Receptor Binding Domain (RBD) in complex with hACE2 (starting point for each simulation was modelled off crystal structures of the SARS-Cov2 Spike-RBD/hACE2 complex),27 which we used to calculate the free energy contributions from electrostatics, polar solvation, van Der Waals interactions, and solvent-accessible surface area (SASA) to infer the free energy of binding between those two proteins.40,41 We quantified the root-mean-squared deviation (RMSD) of the portion of the RBD closest to the hACE2 receptor (residues 397 to 512) for each of our replicates to confirm complex stability (Supplementary Figure 2). Contrary to our expectations, the free energy of binding between the Spike-RBD and the hACE2 appears to have decreased between N0 and N1. In fact, each of the 4 changes either reduced or did not significantly change the free energy of binding (Figure 2A) (this is true for both alternate reconstructions of position 498 in N0). While this was somewhat surprising, it corresponds with recently released in vitro binding measurements remarkably.42 In particular, we see that both alternative reconstructed states for position 498 in N0 clearly improve hACE2 binding affinity in both our simulations and in in vitro functional measurements.42
While overall it is clear these four historical changes reduced binding affinity to hACE2, they did not do so equally: t346R and h/y498Q showed the largest decreases in affinity (Figure 2B). These results demonstrate that, contrary to expectations, evolutionary changes since 2013 did not improve the SARS-Cov2 Spike-RBD’s binding with hACE2. To our knowledge, this is the first demonstration that the SARS-CoV-2’s common ancestor with the RaTG13 lineage may have been capable of binding to hACE2. This has important implications for our understanding of how SARS-CoV-2 evolved to infect humans. First, it suggests that the binding affinity between the Spike-RBD and hACE2 may not be a critical driver of the high degree of infectivity that has been observed during its recent outbreak. Instead, it suggests that tight hACE2 binding may be a latent property of this virus, and that high infectivity may instead have emerged via a distinct set of molecular changes in the SARS-Cov2 genome. Second, it calls into question the presumption of a recent zoonotic origin for this disease; while other molecular components of the current SARS-Cov2 virus may have acquired recent evolutionary changes that promoted its infectivity in humans, it appears that the high affinity for hACE2 was not among them.
If this is the case – that this viral lineage possessed the ability to bind hACE2 with high affinity for at least the past 7 years (Figure 1B) – then why did it not emerge as a public health issue until recently? One possibility is that binding hACE2 by the Spike-RBD is not sufficient, on its own, to infect humans, and that other molecular components first needed to acquire new functions to do so. A second possibility is that this virus may have been capable of infecting human cells for a longer period of time in the past, but that its ancestral form either presented with far fewer symptoms (making it less disruptive and/or noticeable to those infected), or that it was far less infectious (thereby impacting only a small number of people), in either case escaping the notice of public health monitoring (Figure 2C). To test this, a broad and concerted effort to sequence the range of coronaviruses across human populations would need to be conducted, in order to test whether a closely related virus may also be circulating.43–45
Naturally, as an in silico study, these results should be interpreted with some caution. Insofar as they can be validated, however, they are largely consistent with direct functional measurements in the lab.42 Ideally, combinatorial libraries could be constructed46,47 and functionally screened48 in order to glean more detailed insights into the molecular mechanisms underlying the recent evolution of this virus.
Predicting the emergence of highly infectious and virulent diseases, while difficult, is vital for human population health.49 To do so, however, we must take steps to understand how pandemic diseases – such as SARS-Cov2 – emerged as they did, and to understand if and when they acquired the novel molecular functions that enabled their infectivity. In this case, it appears that the SARS-Cov2 Spike-RBD did not recently evolve binding affinity to a human-specific protein. Instead, that function appears to have been latent, making it clear that the evolution of this disease – along with so many other aspects of its etiology – is more complex than expected.
Confirmation of SARS-CoV-2 genome evolution
A phylogenetic analysis of 26 viral genomes was performed to confirm known SARS_CoV_2 ancestors. 24 known enzootic and endemic viruses and the SARS-CoV-2 reference genome and the Pangolin-CoV genome were downloaded from the National Center for Biotechnology Information (NCBI)50 and Lam et al.21 respectfully. Selected sequences were aligned using the Multiple Alignment using Fast Fourier Transform Version 7 (MAFFT) FFT-NS-2 algorithm. 51,52 MAFFT default parameters were used in our alignment, meaning gap penalties were assigned a value of 1.53. PhyML 3.0 was employed to construct a phylogeny of aligned genomes. 53,54 Bayes values ≥ 0.90 were considered statistically significant. The output tree was visualized using the online tool, Interactive Tree of Life (iTOL), and statistically significant clades were examined to validate current knowledge surrounding SARS-CoV-2 evolution.55
We isolated the RbRp domain from the SARS-CoV-2 reference genome using the Basic Local Alignment Search Tool (BLAST).56 First, we performed a tblastn search of the SARS-CoV-2 RbRp reference domain published on NCBI using a BLOSUM 62 matrix, a gap opening penalty of 11 and a gap extension penalty of 1.50 We then used the output of this query to find the specific location of the RbRp domain in the nucleotide SARS-CoV-2 reference genome and isolated this portion of the genome for our analysis. Next, we employed BLASTn t
o isolate the RbRp domain from the pangolin coronavirus and the RaTG13 coronavirus. In this alignment, we used the gap opening penalty of 0 and gap extension penalty of 2.5. We also downloaded the 127 RbRp bat coronavirus sequences published by Joffrin et al.15 Isolated RbRp sequences were aligned using the MAFFT G-INS-I algorithm.51,52 Gap opening and extension penalties were the same as previously described. Finally, we created a phylogeny of the aligned sequences using PhyML 3.0. 53,54 The tree was visualized using iTOL and statistically significant clades were examined.55
Identification of geographic differences in SARS-CoV-2 sequences
To identify geographic and time dependent differences in the SARS-CoV-2 viral genome, we performed a phylogenetic analysis of 479 SARS-CoV-2 sequences obtained from GISAID (Supplementary Table 3).57,58 We arbitrarily selected one sequence per day per country from December 30, 2019 to March 25, 2020. We also included the RaTG13 reference genome and the Pangolin-CoV genome published by Lam et al. in our analysis. 21,50 Selected sequences were aligned, as previously described, using the MAFFT FFT-NS-2 alignment algorithm.51,52 A consensus sequence was constructed from the 479 aligned SARS-CoV-2 sequences using the online tool MSAViewer.59 We validated our consensus sequence by cross referencing the sequence to the SARS-CoV-2 reference genome published by NCBI.50 Following, validation we uploaded the consensus sequence to the online cloud-based informatics platform, Benchling, and extracted the ORFs (<a href=”http://benchling.com” rel=”nofollow”>benchling.com</a>).
To determine common SNPs within the 479 SARS-CoV-2 sequences, the consensus ORFs were aligned using Nucleotide BLAST.56 By default, our analysis used a gap opening penalty of 5, a gap extension penalty of 2 and a mismatch penalty of −3. SNPs were extracted from the BLAST output using the Java module BlastNToSnp (http://dx.doi.org/10.6084/m9.figshare.1425030). We then selected for SNPs present in >5% of the sequences and inputted the SNPs into Benchling to determine whether the SNPs caused silent or missense mutations in the ORFs (<a href=”http://benchling.com” rel=”nofollow”>benchling.com</a>). PhyML 3.0 was employed to construct a phylogenetic tree of the 479 aligned SARS-CoV-2 genomes, the RaTG13 genome and the Pangolin-CoV genome, 53,54 which was visualized using iTOL.55
Ancestral sequence reconstruction of spike glycoprotein receptor binding domain
nBLASTx, run using a BLOSUM 62 matrix, a gap opening penalty of 11 and a gap extension penalty of 1, was employed to extract the Spike glycoprotein from the 479 SARS-CoV-2 sequences and the Pangolin-CoV genome.56 Additional, Spike sequences, including the RaTG13 Spike protein, were obtained directly from NCBI.50 Protein sequences were initially aligned using the Multiple Sequence Alignment by Log-Expectation (MUSCLE) program.60 The optimal parameters for phylogenetic reconstruction analysis were taken from the best-fit evolutionary model selected using the Akaike Information Criterion (AIC) implemented in the PROTTEST3 software,61 and were inferred to be the Jones-Taylor-Thornton (JTT) model62 with gamma-distributed among-site rate variation and empirical state frequencies. Phylogeny was inferred from these alignments using the RaXML v8.2.9 software63 and results were visualized using FigTree v1.4.4 (https://github.com/rambaut/figtree/releases). Ancestral sequence reconstruction was performed with the FastML software64 and further validated independently using the Graphical Representation of Ancestral Sequence Predictions (GRASP) software.65 Statistical confidence in each position’s reconstructed state for each ancestor determined from posterior probability; any reconstructed positions with less than 95% posterior probability was considered ambiguous, and alternate states were also tested.
Mutagenesis of ancestral proteins
To understand the evolutionary importance of sequence changes observed between ancestral, zoonotic, and SARS-CoV-2 spike protein sequences, in silico mutagenesis and binding energy studies were performed. A previously constructed x-ray crystallography structure for the complex between the receptor binding domain (RBD) of the SARS-CoV-2 spike protein and the human hACE2 receptor were obtained from RCSB (accession number 6M0J). Utilizing PyMOL mut
agenesis wizard, 66 the four missense mutations (R346t, A372t, Q498h or Q498y, H519n) identified between the N0 and N1 sequences were introduced into the SARS-CoV-2 RBD sequence, replicating the sequence of the putative ancestral zoonotic (N0) sequence. In addition to the N1 and N0 structures, additional structures were developed in a similar fashion, selectively including each of the 4 mutations to represent all of the possible combinations that these mutations may have existed throughout evolutionary time
Molecular dynamics simulation of Spike-RBD-hACE2 interactions
Molecular interactions were characterized with molecular dynamics simulations using Gromacs, TIP3P waters and CHARM07 force-field parameters for proteins. For each condition, three replicate 10 ns simulations were run, starting from crystal structures or structural models. Historical mutations were introduced and energy-minimized before MD simulation. Each system was solvated in a cubic box with a 10 Å margin, then neutralized and brought to 150 mM ionic strength with sodium and chloride ions. This was followed by energy minimization to remove clashes, assignment of initial velocities from a Maxwell distribution, and 1 ns of solvent equilibration in which the positions of heavy protein and DNA atoms were restrained. Production runs were 50 ns, with the initial 10 ns excluded as burn-in. The trajectory time step was 2 fs, and final analyses were performed on frames taken every 12.5 ps. We used TIP3P waters and the CHARM07 FF03 parameters for proteins, as implemented in GROMACS 22.214.171.124 Analyses were performed using VMD 126.96.36.199 GROMACS output was uploaded into Visual Molecular Dynamics (VMD) for Root-Mean Squared Deviation (RMSD) Analysis using the RMSD trajectory tool (ref). After discovering large deviations in RMSD values for the full RBD, which we attributed to noise at the ends of the RBD, we isolated our analysis to residues 397 to 512 of the RBD.
Measurement of binding energies
Next, we measured the binding energies between residues 397 to 512 and the ACE2 receptor using g_mmpbsa, a program which employs Molecular mechanics Poisson–Boltzmann surface area (MMPBSA) calculations to determine binding energy. Van der Waal forces, polar solvation energy, apolar solvation energy and SASA energy were calculated every 0.25 ns using a gridspace of 0.5 and a solute dielectric constant of 2. The output of the three replicates was amalgamated and binding energy was calculated using the bootstrap analysis (n = 2000 bootstraps) published by Kumari et al.40,41 We then characterize the genetic effect of each mutation (on average) and assessed whether there were any statistically significant epistatic interactions using established methods.46,47
↵+ Co-first authors
Michael Novakhov – SharedNewsLinks℠