key: cord-310477-vniokol0
authors: Pontes, Camila; Ruiz-Serra, Victoria; Lepore, Rosalba; Valencia, Alfonso
title: Unraveling the molecular basis of host cell receptor usage in SARS-CoV-2 and other human pathogenic β-CoVs
date: 2020-08-21
journal: bioRxiv
DOI: 10.1101/2020.08.21.260745
sha: 
doc_id: 310477
cord_uid: vniokol0

The recent emergence of the novel SARS-CoV-2 in China and its rapid spread in the human population has led to a public health crisis worldwide. Like in SARS-CoV, horseshoe bats currently represent the most likely candidate animal source for SARS-CoV-2. Yet, the specific mechanisms of cross-species transmission and adaptation to the human host remain unknown. Here we show that the unsupervised analysis of conservation patterns across the β-CoV spike protein family, using sequence information alone, can provide rich information on the molecular basis of the specificity of β-CoVs to different host cell receptors. More precisely, our results indicate that host cell receptor usage is encoded in the amino acid sequences of different CoV spike proteins in the form of a set of specificity determining positions (SDPs). Furthermore, by integrating structural data, in silico mutagenesis and coevolution analysis we could elucidate the role of SDPs in mediating ACE2 binding across the Sarbecovirus lineage, either by engaging the receptor through direct intermolecular interactions or by affecting the local environment of the receptor binding motif. Finally, by the analysis of coevolving mutations across a paired MSA we were able to identify key intermolecular contacts occurring at the spike-ACE2 interface. These results show that effective mining of the evolutionary records held in the sequence of the spike protein family can help tracing the molecular mechanisms behind the evolution and host-receptors adaptation of circulating and future novel β-CoVs. Significance Unraveling the molecular basis for host cell receptor usage among β-CoVs is crucial to our understanding of cross-species transmission, adaptation and for molecular-guided epidemiological monitoring of potential outbreaks. In the present study, we survey the sequence conservation patterns of the β-CoV spike protein family to identify the evolutionary constraints shaping the functional specificity of the protein across the β-CoV lineage. We show that the unsupervised analysis of statistical patterns in a MSA of the spike protein family can help tracing the amino acid space encoding the specificity of β-CoVs to their cognate host cell receptors. We argue that the results obtained in this work can provide a framework for monitoring the evolution of SARS-CoV-2 specificity to the hACE2 receptor, as the virus continues spreading in the human population and differential virulence starts to arise.

135 spike protein sequences which were manually reviewed to annotate information on viral strain, host organism and cell receptors using information extracted from the NCBI Protein database, UniprotKB and the literature. Detection of SDPs and spike protein subfamilies. The S3Det method (15) uses multiple correspondence analysis (MCA) to identify differentially conserved positions and sequence subfamilies within a given MSA. MSA positions that follow the subfamily segregation are defined as specificity determining positions (SDPs) of the family. Here, the unsupervised mode of S3Det was used in a two-level decomposition analysis to identify SDPs linked to the spike protein family segregation between and within β-CoV subgroups. Phylogenetic analysis was performed both on the full β-CoV MSA and for individual subfamilies using the PhyML method (16) with default parameters. Coevolution analysis. Coevolving MSA positions were identified by computing the MI-APC (17) . MI-APC is a mutual information (MI)-based score corrected by the average product correction (APC) of the background noise and phylogenetic signal. To ensure robust statistics, MSA columns were filtered according to percentage of gaps (<= 50%) and Shannon entropy (>= avg -std), computed as follows: where x is the number of SDPs observed in a given domain of interest, k is the number of SDPs in the test set, m is the length of the domain of interest and n + m is the size of the protein. The analyses were performed both at the subunit and domain level. Annotations on subunits and domain boundaries of the SARS-CoV-2 spike glycoprotein were retrieved from the literature (18) and mapped to the rest of human β-CoVs spike sequences used in this analysis (SARS-CoV, MERS, OC43 and HKU1) based on the MSA. Evolutionary rate analysis. Per-site evolutionary rates (EV) were computed by Rate4Site (19) using as input the full MSA of the β-CoV spike family and phylogenetic tree. SARS-CoV-2 sequence was set as the reference sequence. Raw rate values were used to compute relative rates by normalizing the values to the mean of 1 (20) . Structural dataset. All the structures employed in this study were retrieved from the PDB (https://www.rcsb.org) using the following PDB codes: 6VSB (7) , 6VXX (21) , 6LZG (22) , 5X5B (23) , 5X58 (23) , 2AJF (24) , 6OHW (25) , 6NZK (25) , 5X5F (23) , 4L72 (26) and 5I08 (27) .

A phylogenetic analysis was performed based on a MSA of the full-lenght spike sequences from SARS-CoV-2 and other representative β-CoVs. The resulting tree, shown in Figure 1A , reflects the taxonomic classification of β-CoVs into five subgenera, namely Sarbecovirus, Hibecovirus, Nobecovirus, Embecovirus and Merbecovirus (29) . A first S3Det analysis was performed on individual protein subfamilies from the Sarbecovirus, Embecovirus and Merbecovirus subgroups. In the case of Sarbecovirus, we identified three different clusters, two of which correspond to known Sarbecovirus clades (30) . Notably, in contrast to what observed based on phylogenetic analysis, both SARS-CoV-2 and RaTG13 are clustered together with members of Sarbecovirus clade 1, which includes human SARS-CoV and bat SARS-like sequences ( Figure 1B-C) . This result is confirmed by the analysis of similarity scoring matrices (Supplementary Figure S1) where SARS-CoV-2 and SARS-CoV sequences cluster together and are closer to the green cluster based on the similarity of SDPs, but not when the full-length sequence is considered. Within the Embecovirus group, we identified three clusters: a first cluster corresponding to murine CoVs (MHV), a second containing rat (RtCoVs) and human CoVs (CoV-HKU1), and a third cluster containing the human CoV hCoV-OC43 and other mammalian CoVs (Supplementary Figure S2 -A,C). Within the Merbecovirus, we identified four clusters: a first cluster containing MERS and MERS-like bat CoVs, a second cluster containing hedgehog and bat CoVs, a third cluster containing HKU5 bat CoVs, and a fourth cluster containing one single bat CoV isolate (KW2E-F93) from the Nycteria species (Supplementary Figure S2-B,D) . The SDPs associated with the subfamily segregation within these three subgenera display a strong domain enrichment within the spike S1 subunit ( Figure 2 , Supplementary Table S1 ). Specifically, the SDPs (n = 28) of Sarbecovirus subfamilies, containing both SARS-CoV and SARS-CoV-2 spike sequences, are enriched in the RBD whereas in hCoV-OC43 and hCoV-HKU1, the human-pathogenic species belonging to Embecovirus, SDPs (n = 12) fall mainly in an upstream region of the spike, with a significant enrichment in the NTD. Also in MERS-CoV we observe a significant enrichment within the S1 subunit. However, as shown in Figure 2 , the SDPs (n = 33) are almost evenly distributed across the NTD and RBD regions, with a significant enrichment in the latter (Supplementary Table S2) . Notably, the distribution of the SDPs shows a clear relationship with the cell receptor usage observed among β-CoVs. In particular, hCoV-HKU1 and hCoV-OC43 are known to bind sialic acid receptors on the host cells via the spike NTD (31) , SARS-CoV and SARS-CoV-2 recognize the ACE2 receptors via the RBD (32, 33) , while MERS-CoV has been reported to use both DPP4 (34) and sialic acid receptors (35) via the RBD and NTD domains, respectively. In summary, the SDPs found within these β-CoV subgenera define a specific region of the receptor binding domains: they are part of, or in direct contact with, the ACE2 interacting surface ( A second S3Det analysis was performed on the full β-CoV MSA. As shown in Figure 1A -D, the identified sequence subfamilies are consistent with the phylogenetic classification of Betacoronavirus into five subgenera (29) . In this case, the SDPs linked to the full β-CoV spike family segregation are mainly located in the S2 protein subunit, with a statistically Figure S4) . Collectively, these results indicate that SDPs capture the functional diversification observed within the individual protein subfamilies, whereby host-cell receptor specificity arises in the context of a structural framework that is specific to each β-CoV phylogenetic group.

Structural and mutagenesis studies have shown that the spike RBD of Sarbecovirus contains all the necessary information for host receptor binding and that a few amino acid substitutions in this region can lead to efficient cross-species transmission (30, 37) . Binding to ACE2 is clade-specific and occurs at the carboxy-terminal region of the RBD, by an extended concave loop subdomain which forms the interaction interface with the ACE2 N-terminal helix. Notably, both the ACE2-contacting residues and the surrounding amino acids, collectively referred to as the receptor binding motif (RBM), are required to impart human receptor usage within the Sarbecovirus lineage (30) . Consistently, we observe that several SDPs associated with the Sarbecovirus subfamily fall within the RBM, i.e., Y451/Y438, L461/L448, T470/N457, C488/C474, P491/P477, G496/G482, G502/G488, P507/P493, Y508/Y494 on SARS-CoV-2/SARS-CoV sequences, respectively ( Figure 3 ). Of these, Y494 has been previously reported as critical for ACE2 binding (38) . Two SDPs, namely G496 and G502, fall within the receptor interface forming two hydrogen bonds with the ACE2 K353. Other SDPs, such as L461, T470, C488, P491 and G496 make direct contact with ACE2 contacting residues. Hence, these SDPs are likely to play an important role in ACE2 binding by affecting the local orientation of ACE2 contacting residues. This hypothesis is further supported by results from in silico mutagenesis and coevolution analysis (details in next section). Specifically, we tested the effect of amino acid mutations across the RBM by mutating the SARS-CoV-2 sequence to the consensus spike sequence from clade 2, i.e. a clade known to be incompatible with hACE2 usage (30) . Notably, while most substitutions are predicted to have a destabilizing effect on the Spike-hACE2 complex, mutations of SDP residues are predicted to have a significantly larger impact compared to non-SDP residues (Figure 4) .

Computational methods exploiting coevolution signals in MSAs of protein families are widely used to infer features such as molecular interactions and functional sites (15, (39) (40) (41) . Such signals arise from the specific adaptation between correlated amino acid sites, where changes in one site are potentially compensated by changes in the other. In the case at hand, coevolution signatures are used as markers for the study of the physical interactions occurring between different sites of the spike as well as between the spike and their cognate host cell receptor ACE2. As it can be observed in Supplementary Figure S6 , the strongest intramolecular coevolution signal, considering the top-500 predictions, is observed over the RBD region of the spike (the overall precision of the method is reported in Supplementary Figure S7 ). Figure 5 shows in detail the RBM region, which presents above average precision and recall values of 46% and 4.6%, respectively. Among the SDPs (highlighted in green), the precision is even higher, around 62%, and the recall is 8.7%. Notably, 18% of SDPs found within the Sarbecovirus subgenus show a coevolution signal with ACE2 contacting residues, namely Y489 (coevolving with C488, P491), Q493 (coevolving with L461, T470, C488, P491) and Q498 (coevolving with G496). These three positions are hubs on the interface with ACE2, making direct contact with positions T27, F28, K31 and Y83, positions K31, H34 and E35, and positions D38, Y41, Q42 and L45 of ACE2, respectively (22) . Particularly, position Q493/N479 in SARS-CoV-2/SARS-CoV has been described to be critical for high affinity binding of both SARS-CoV and SARS-CoV-2 to ACE2 (42, 43) . Furthermore, it is interesting to notice that the C488/C474 SDP in SARS-CoV-2/SARS-CoV is an important position for the stability of the RBM as a whole and a complete loss of hACE2 binding in vitro has been described when this position is mutated to Alanine in SARS-CoV (44) .

We next performed a coarse-grained coevolutionary analysis on a concatenated MSA containing eight spike proteins and their cognate ACE2 receptors. Contact predictions were obtained by computing the MI-APC score for every inter-protein pair of alignment positions, considering the RBD region of the spike protein and the whole ACE2. Interestingly, three RBM positions were found coupled to ACE2 positions among the top-10 MI-APC scores (Supplementary Figure S8) . Specifically, the ACE2 residue H34 was coupled to L455, S494 and Q498 on the spike protein. Additionally, the spike positions R346 and L455 were coupled to ACE2 residues Q24, N61, Q81, M82 and to E23, T27 and H34, respectively. Among these predictions, T27-L455, H34-L455 and H34-S494 correspond to true contacts within 8A distance, while the couplings between R346 and ACE2 positions could be related to long-range effects. Notably, positions E23 and H34 have been described as crucial to SARS-CoV binding to ACE2 (45) . Also, L455 was described as important for the stabilisation of a binding hotspot between SARS-CoV-2 and ACE2 (43) . Interestingly, a recent study reported that ACE2 variants E23K and T27A are more susceptible to SARS-CoV-2, while variant H34R decreases SARS-CoV-2 affinity (46) . These results, despite the limited number of sequences in the concatenated alignment, point to at least three specific RBM viral positions (L455, S494 and Q498 in SARS-CoV-2) likely adapted to their species-specific counterparts.

Since the beginning of the COVID-19 outbreak and the isolation of the SARS-CoV-2 virus, laboratories around the world are continuously isolating viral genomic sequences with unprecedented speed, enabling nearly real-time data sharing of more than 80,000 genomic sequences so far (47) . After discarding partial sequences, a multiple sequence alignment was built based on a total of 44,812 SARS-CoV-2 spike sequences isolated from human samples in 96 countries. Our analysis of missense amino acid variations confirmed earlier reports (48) that most mutations occur within the S1 subunit, with a dominant variant observed at position 614, where more than 75% of samples carry the D614G mutation, followed by mutations L4F and R21I appearing in 337 and 307 samples, respectively. Within the RBD, the top frequent mutations are S477N, T478I and P479S, found in 68, 63, and 53 samples, respectively. Notably, these positions fall within the RBM, forming a surface-exposed loop that is proximal to the ACE2 binding surface, and that is absent in Sarbecovirus clades 2 and 3 (30) . Although previous experiments indicate that this loop is per se not sufficient to impart ACE2 receptor usage, deletions of this region are associated with reduced spike expression and loss of cell entry (30) . Hence, mutations in this region are expected to impact the stability of the protein, rather than its affinity to the receptor. Finally, perhaps not surprisingly, the frequency of variants at SDPs is very low, with an overall variability that is comparable to that observed in ACE2-contacting residues (Supplementary Figure S9) . Mutations are observed in 7 out 14 SDPs. With the exception of T385A and Y508, that are observed in 2 and 4 samples, respectively, all other mutations are only present in one sample. Mostly, these variants are predicted to be neutral or to destabilize the binding to the receptor (Supplementary Figure S10B) . In summary, the low frequency of mutations in SDPs fulfils the requirement to preserve a functional role of those positions, providing an additional evidence of their involvement in the interaction with the host cell receptor.

The relationship between protein family segregation and their functional organisation has been extensively investigated for decades and a variety of computational methods have been developed to infer their evolutionary link at the residue level (49) (50) (51) . It is therefore relatively straightforward to identify the amino acid positions that modulate the functional specificity of a given enzyme towards a substrate or cofactor (15) or the binding specificity of a protein-ligand or protein-protein interaction (52, 53) by the analysis of the differential conservation patterns within the MSA of a protein family. Here we apply this concept to the analysis of the SARS-CoV-2 spike in the context of a MSA of homologous sequences belonging to the β-CoV genus. The SDP analysis is based on a vectorial representation of protein sequences and amino acid positions in a multidimensional space to simultaneously identify the family segregation and the residue positions that better explain the sources of variation of the family (15) . On one hand, the analysis performed at the β-CoV family level led to the identification of five sequence subfamilies, reflecting the known phylogenetic classification of Betacoronavirus into five subgenera (29) ( Figure 1A ). On the other hand, the analysis performed on individual β-CoV subgenera, i.e. Sarbecovirus, Merbecovirus and Embecovirus subgroups, allowed a fine-grained classification into subfamily clusters that clearly reflect the functional diversification of the spike protein family, that is, the specificity to different host-cell receptors (Figure 2-3) . Indeed, both the clustering and domain enrichment results of the associated SDPs consistently reproduce the known cell receptor specificities observed across the different β-CoV lineages (26, (54) (55) (56) . At the level of the Sarbecovirus group, for example, both SARS-CoV-2 and RaTG13 are clustered together with SARS-CoV and other SARS-like sequences from bats ( Figure 1C) , reflecting the ability of these members of the Sarbecovirus group to bind the ACE2 cell receptor (57) . Notably, the proximity of SARS-CoV-2 and RaTG13 and other SARS-CoV sequences based on key SDP positions is different to what seen based on full sequence phylogenetic analysis, where they form distinct clades. A proximity that is driven by their shared SDPs and it is interpreted here in terms of their shared ability to bind the human ACE2 receptor. As it is often the case, functional constraints arise from the requirement of maintaining the interaction of proteins with other macromolecules or ligands. Such constraints translate into specific roles and properties of individual amino acids, or protein sites. In the case at hand, the analysis of the physicochemical, structural and conservation properties of the SDPs of the different subfamilies highlights a pattern that is typical of protein functional sites, as they show high conservation across the protein family, are solvent exposed, and are enriched in the receptor binding domains (Supplementary Figures S3-5 and Supplementary Table S2) . Hence, in order to assess the role of the SDPs in mediating ACE2 receptor usage across the Sarbecovirus group, we set up an in silico mutagenesis study and analysed the effect of amino acid mutations across the RBM and their impact on ACE2 binding. Notably, while our results are in line with previous observations (30) , they point to specific positions across the RBM that might exert a critical role, either by engaging the receptor through direct intermolecular contacts or by affecting the local orientation of ACE2 contacting residues and the stability of the RBM as a whole. Collectively, these results point to a key role of SDPs in mediating host cell receptor specificity across β-CoVs and provide, at the same time, a framework for monitoring the evolution of the SARS-CoV-2 specificity to hACE2, as well as the emergence of novel potential cross-species transmission events. As such, it is important to notice that from the analysis of amino acid variations across the circulating SARS-CoV-2 virus, SDPs tend to mutate with a very low frequency, similar to what is seen at ACE2 contacting sites (Supplementary Figure S9) (58) . This is of relevance, as our results suggest that mutations in SDPs can significantly impact the receptor-binding ability of the spike. Furt hermore, the experience in other scenarios has shown that mutating SDPs is, in general, sufficient to transform the properties between two groups of proteins of the same family, i.e. the interchange of the residues occupying the SDPs between two families implies a change in the associated biological properties (59) (60) (61) (62) (63) . Notable examples include the production of switch-of-function mutants of small GTPases with changed selectivity, or the change of transport specificity between MIP channel proteins by few amino acid substitutions (60, 62) . In line with this reasoning, it can be argued that other members of the Sarbecovirus group might have the potential to acquire ACE2 binding ability, as they share substantial similarity in terms of SPDs. This is especially the case of members of the Sarbecovirus green cluster, which despite being phylogenetically distant to SARS-CoV-2, display identical residues in 10 out of 14 SDP positions within the RBD, making them potential candidates for new human infections.

In conclusion, the results presented here show that the identification of evolutionary patterns based on the analysis of sequence information alone can provide meaningful insights on the molecular basis of host-pathogen interactions and adaptation. We believe that both the methodology and results presented in this work can provide the basis for follow-up studies analysing the potential routes of mutations that could lead to new adaptation to human hosts and ultimately contribute to better understanding and monitoring of events that are critical to public health concerns worldwide.

The SARS epidemic in Hong Kong: what lessons have we learned?

Outbreak of Middle East respiratory syndrome coronavirus in Saudi Arabia: a retrospective study

An interactive web-based dashboard to track COVID-19 in real time

Genetic Recombination, and Pathogenesis of Coronaviruses

Fusion mechanism of 2019-nCoV and fusion inhibitors targeting HR1 domain in spike protein

Evidence for a common evolutionary origin of coronavirus spike protein receptor-binding subunits

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor

The proximal origin of SARS-CoV-2

Structure, Function, and Evolution of Coronavirus Spike Proteins

Evolutionary origins of the SARS-CoV-2 sarbecovirus lineage responsible for the COVID-19 pandemic

Basic local alignment search tool

CD-HIT: accelerated for clustering the next-generation sequencing data

MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization

Protein interactions and ligand binding: from protein subfamilies to functional specificity

New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0

Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction

Inhibition of SARS-CoV-2 (previously 2019-nCoV) infection by a highly potent pan-coronavirus fusion inhibitor targeting its spike protein that harbors a high capacity to mediate membrane fusion

Rate4Site: an algorithmic tool for the identification of functional regions in proteins by surface mapping of evolutionary determinants within their homologues

Measuring evolutionary rates of proteins in a structural context

Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein

Structural and Functional Basis of SARS-CoV-2 Entry by Using Human ACE2

Cryo-EM structures of MERS-CoV and SARS-CoV spike glycoproteins reveal the dynamic receptor binding domains

Structure of SARS coronavirus spike receptor-binding domain complexed with receptor

Structural basis for human coronavirus attachment to sialic acid receptors

Structure of MERS-CoV spike receptor-binding domain complexed with human receptor DPP4

Pre-fusion structure of a human coronavirus spike protein

Estimating maximum likelihood phylogenies with PhyML

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

Functional assessment of cell entry and receptor usage for SARS-CoV-2 and other lineage B betacoronaviruses

Human coronaviruses OC43 and HKU1 bind to 9--acetylated sialic acids via a conserved receptor-binding site in spike protein domain A

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Angiotensin-converting enzyme 2 is a functional receptor for the SARS coronavirus

Identification of residues on human receptor DPP4 critical for MERS-CoV binding and entry

Identification of sialic acid-binding function for the Middle East respiratory syndrome coronavirus spike glycoprotein

Slow protein evolutionary rates are dictated by surface-core association

A SARS-like cluster of circulating bat coronaviruses shows potential for human emergence

The SARS Coronavirus S Glycoprotein Receptor Binding Domain: Fine Mapping and Functional Characterization

Sequence co-evolution gives 3D contacts and structures of protein complexes

Emerging methods in protein co-evolution

Conservation of coevolving protein interfaces bridges prokaryote-eukaryote homologies in the twilight zone

Receptor and viral determinants of SARS-coronavirus adaptation to human ACE2

Structural basis of receptor recognition by SARS-CoV-2

A 193-amino acid fragment of the SARS coronavirus S protein efficiently binds angiotensin-converting enzyme 2

Identification of critical determinants on ACE2 for SARS-CoV entry and development of a potent entry inhibitor

Human ACE2 receptor polymorphisms predict SARS-CoV-2 susceptibility https

Buckland-Merrett, Data, disease and diplomacy: GISAID's innovative contribution to global health

Characterizations of SARS-CoV-2 mutational profile, spike protein stability and viral transmission

Evolution of function in protein superfamilies, from a structural perspective

A method to predict functional residues in proteins

Protein function prediction using local 3D templates

An evolutionary trace method defines binding surfaces common to protein families

Automatic methods for predicting functionally important residues

Molecular basis of binding between novel human coronavirus MERS-CoV and its receptor CD26

Bat origins of MERS-CoV supported by bat coronavirus HKU4 usage of human receptor CD26

Receptor usage and cell entry of bat coronavirus HKU4 provide insight into bat-to-human transmission of MERS coronavirus

Mutations from bat ACE2 orthologs markedly enhance ACE2-Fc neutralization of SARS-CoV-2 https

Predicting functionally important residues from sequence conservation

Evolution of protein kinase substrate recognition at the active site

Switch from an Aquaporin to a Glycerol Channel by Two Amino Acids Substitution

Effector recognition by the small GTP-binding proteins Ras and Ral

Switch-of-Function Mutants Based on Morphology Classification of Ras Superfamily Small GTPases

Evolution-guided discovery and recoding of allosteric pathway specificity determinants in psychoactive bioamine receptors

The authors are grateful to all other members of the Computational Biology group for their valuable feedbacks and useful discussions.

This work has received funding from the EXSCALATE4CoV project, from the European Union's Horizon 2020 Research and Innovation Programme, under grant agreement N. 101003551.

The authors have declared no competing interest.