Quantifying the immunological distinctiveness of emerging SARS-CoV-2 variants in the context of prior regional herd exposure

Abstract The COVID-19 pandemic has seen the persistent emergence of immune-evasive SARS-CoV-2 variants under the selection pressure of natural and vaccination-acquired immunity. However, it is currently challenging to quantify how immunologically distinct a new variant is compared to all the prior variants to which a population has been exposed. Here, we define “Distinctiveness” of SARS-CoV-2 sequences based on a proteome-wide comparison with all prior sequences from the same geographical region. We observe a correlation between Distinctiveness relative to contemporary sequences and future change in prevalence of a newly circulating lineage (Pearson r = 0.75), suggesting that the Distinctiveness of emergent SARS-CoV-2 lineages is associated with their epidemiological fitness. We further show that the average Distinctiveness of sequences belonging to a lineage, relative to the Distinctiveness of other sequences that occur at the same place and time (n = 944 location/time data points), is predictive of future increases in prevalence (Area Under the Curve, AUC = 0.88 [95% confidence interval 0.86 to 0.90]). By assessing the Delta variant in India versus Brazil, we show that the same lineage can have different Distinctiveness-contributing positions in different geographical regions depending on the other variants that previously circulated in those regions. Finally, we find that positions that constitute epitopes contribute disproportionately (20-fold higher than the average position) to Distinctiveness. Overall, this study suggests that real-time assessment of new SARS-CoV-2 variants in the context of prior regional herd exposure via Distinctiveness can augment genomic surveillance efforts.


Introduction
To date, over 10 billion COVID-19 vaccine doses have been administered globally (1), with over 200 million individuals fully vaccinated in the United States (2). Recent studies have confirmed that natural immunity (i.e. immunity gained through prior infection) is also highly protective and may even provide more durable protection than vaccination alone (3)(4)(5)(6)(7)(8)(9)(10)(11)(12). Given that over 400 million COVID-19 cases have been reported worldwide (with over 78 million cases in the United States) (1), it is likely that both vaccination-acquired immunity and natural immunity play important roles in the evolution of new SARS-CoV-2 variants.
Throughout the course of the COVID-19 pandemic, SARS-CoV-2 has evolved to generate new variants which harbor unique constellations of mutations (substitutions, deletions, and insertions). Some of these variants are designated as Variants of Concern (VOCs) based on evidence for increased transmissibility, increased disease severity, or reduced neutralization by vaccine-elicited sera or authorized monoclonal antibody treatments. Such variants include Alpha (B.1.1.7 and Q lineages per PANGO classification), Beta (B.1.351 and descendants), Gamma (P.1 and descendants), Delta (B.1.617.2 and AY lineages), and most recently Omicron (B.1.1.529 and BA lineages) (13). As new SARS-CoV-2 variants evolve, it is im- portant to estimate their likelihood of evading existing regional herd immunity and potentially transmitting highly at the community level. While the main evidence that a variant evades immune response typically relies on laboratory assays and epidemiological evidence (14)(15)(16)(17)(18)(19)(20), complementary approaches that provide initial evidence, as soon as new sequences are reported, could enable earlier response.
Here, we define a new metric "Distinctiveness" to capture the proteome-level novelty of emerging SARS-CoV-2 sequences against all the documented regional lineages. Distinctiveness aims to quantify previous herd exposure to viral sequences that are similar to the current sequence, capturing an important factor of viral epidemiological fitness. This approach views viral evolution through a new lens that considers the pressure to evolve new strains harboring protein content to which communities have not previously been exposed. We show that the same lineage can have different Distinctiveness values simultaneously in different countries, as well as different Distinctiveness-contributing positions. We find that the relative Distinctiveness of emergent SARS-CoV-2 lineages is associated with their epidemiological fitness, as defined by the change in the lineage prevalence. We also show that epitope positions contribute disproportionately to Distinctiveness.

"Distinctiveness" as a metric to capture novelty of emerging SARS-CoV-2 sequences
Understanding the immunological novelty of a SARS-CoV-2 strain for a given population needs to take into consideration which sequences were previously seen at a regional level and for which there might exist population-level immunity. Here, we introduce a new metric "Distinctiveness" of a given SARS-CoV-2 sequence based on comparison against all available sequences previously collected from the same region. Specifically, Distinctiveness is defined as the average distance at the amino-acid level between a sequence and all prior sequences (Figure 1; see the "Methods" section). Distinctiveness can be computed at the global level or at a regional level for any chosen time period. Below, we compare Distinctiveness of the VOCs with contemporary sequences and investigate the relationship between Distinctiveness of a sequence and the change in its regional prevalence. For comparison, we also report the "Mutational load" of the same sequences. Mutational load is simply defined as the number of mutations in the new sequence compared with the ancestral reference sequence (GenBank: MN908947.3), and as such it does not account for the entirety of SARS-CoV-2 evolution or the local prevalence of sequences.
We computed mutational load and Distinctiveness during the emergence of the VOCs in the country of their emergence. Both mutational load and Distinctiveness values of VOC sequences were significantly higher than contemporary lineages (Supplementary Figures S1 and 2). For example, we consider the emergence of the Delta variant in India during January 2021. Both mutational load and Distinctiveness of the Delta variant in India were significantly higher than that of the other contemporary lineages (Figure 2a). This raises the question of whether Delta variant sequences were also competitive in other countries. We considered the example of Brazil, where the Gamma variant was dominant prior to the arrival of the Delta variant ( Figure 2b). Whereas the mutational load of the Delta variant was comparable to those of contemporary lineages, the Distinctiveness of the Delta variant  This position has not contributed to the Delta variant's Distinctiveness as it has been highly prevalent globally (i.e. present in over 99% of SARS-CoV-2 genomes deposited in GISAID) since June 2020 (15,21,22). Brazil, on the other hand, experienced a large wave of cases dominated by the Gamma variant before the arrival of the Delta variant. Here, in addition to the same 10 Spike protein mutations that were observed in India ( Supplementary Figure S3a), there were 11 other positions that further contributed to its regional Distinctiveness (Supplementary Figure S3b). These additional positions correspond to known Gamma lineage-defining mutations (L18F, T20N, P26S, D138Y, R190S, K417T, E484K, N501Y, H655Y, T1027I, and V1176F).

Relative Distinctiveness of emergent SARS-CoV-2 lineages is associated with their epidemiological fitness
In order to examine a possible relationship between Distinctiveness and epidemiological fitness of SARS-CoV-2 lineages, we assessed the correlation between Distinctiveness and change in prevalence for all circulating lineages (grouped as the VOCs and a single group combining all non-VOCs) in 78 geographical regions (27 countries and 51 US states). We find that the relative Distinctiveness of emergent SARS-CoV-2 lineages is associated with their change in lineage prevalence over eight weeks (Figure 3a and Supplementary Figure S4a) (Pearson r = 0.75). In comparison, mutational load was found to have a lower association with change in prevalence (Pearson r = 0.53). We further find that the average Distinctiveness of a lineage in a country/time window can predict future increases in prevalence (Figure 3a, AUC = 0.88 [95% confidence intervals (CIs) 0.86 to 0.90], for predicting a greater than 20 percentage point increase in local prevalence; Supplementary Figure S3b).

Positions that constitute epitopes contribute disproportionately to Distinctiveness compared to the overall SARS-CoV-2 proteome
Since Distinctiveness is intended to capture the fitness of a sequence in the context of previous herd exposure to similar sequences, we next investigated Distinctiveness in the context of known immunogenic positions. Specifically, we analyzed the Distinctiveness of only Spike protein positions and for 66 epitope positions, previously associated with neutralizing antibody binding or therapeutic agent binding (Figure 3b). We found that Distinctiveness determined using only the 66 epitope positions was still correlated with future changes in lineage prevalence (Pearson r = 0.67). Additionally, we found that both Spike positions (average Distinctiveness of 0.007/position) and the 66 epitope positions (average Distinctiveness of 0.061/position) exhibit increased contributions to overall Distinctiveness (average Distinctiveness of 0.003/position).

Discussion
The Distinctiveness metric defined here provides an intuitive quantification of the extent to which any viral sequence differs from other sequences that circulated previously, within the same geographical region. As such, it captures both the emergence of new amino-acid substitutions (e.g. D614G) (23) and deletion of sequence regions that may be involved in antibody recognition (16,24), both of which can affect viral sequence immunogenicity and infectivity. We find that Distinctiveness can predict future changes in local prevalence of newly circulating lineages, suggesting that Distinctiveness could contribute to accurate and early identification of newly circulating lineages that are likely to outcompete other contemporary lineages. For example, analyzing diversity in the Distinctiveness at the US statelevel for the Omicron variant, there are high Distinctiveness sequences in Idaho (Figure 4), warranting a future investigation of sub-regional Distinctiveness within variants and their determinants.
Host immunity against SARS-CoV-2 is largely derived from two sources: vaccination and prior infection. All authorized COVID-19 vaccines utilize the Spike protein sequence from the ancestral Wuhan strain, with a slight modification (substitution of two prolines at positions 986 to 987) to stabilize the prefusion state of the protein product. These vaccines have demonstrated high effectiveness in clinical trials and various real-world studies (17,(25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40), including against most VOCs with the notable recent exception of reduced effectiveness against the Omicron variant (41,42). With over 10 billion vaccine doses administered around the world, it is likely that vaccination-elicited immunity (i.e. antibody and T cell responses against the ancestral Spike protein sequence) acts as a considerable evolutionary pressure on SARS-CoV-2 (1). The importance of natural immunity as an evolutionary pressure is highlighted by several recent studies demonstrating that prior infection confers robust and durable protection against future infection (3)(4)(5)(6)(7)(8)(9)(10)(11)(12). Furthermore, the approach described here can be readily extended to include a correction for the durability of immunity, for example, by reduced contributions to the Distinctiveness calculation of sequences based on their collection date. We suggest that any newly emerging lineage with a combination of sequence modifications that distinguish it from the ancestral strain and VOCs that have circulated widely (or at high prevalence in a given geographic region) should be monitored closely for their potential to drive future surges.
This study has a few limitations. First, we emphasize that the Distinctiveness metric is intended as an intuitive initial evaluation of the novelty of sampled SARS-CoV-2 sequences. By design, it provides a quantification of prior herd exposure, which is a key contributor to population level immunity. However, there are additional factors, such as the functional implications of mutations and stochasticity, that determine whether a new lineage will spread widely that are not captured by our metric and that can only be captured by subsequent lab assays or epidemiological studies. Future work could combine Distinctiveness with such additional contributing factors to develop a more robust predictor of lineage epidemiological fitness. Second, SARS-CoV-2 genomic epidemiology is unfortunately impacted by major geographic and temporal sequencing biases. Over 55% of SARS-CoV-2 genome sequences in GISAID were isolated from infected patients in the United States or the United Kingdom, and the number of cases subjected to whole genome sequencing increased massively starting at the end of 2020. Undersampling of SARS-CoV-2 genomes in other regions and/or during earlier months of the pandemic could impact our estimations of lineage Distinctiveness. Future analysis will include SARS-CoV-2 genomes from complementary databases such as the National Center for Biotechnology Information (43). Third, although we suggest a cut-off for Distinctiveness that can be used to monitor future emerging lineages (Figure 3), it is not clear that this cut-off will remain optimal. The future of SARS-COV-2 evolution is uncertain, and may involve smaller changes to the sequence that necessitate a lower cut-off in Distinctiveness, or a more sensitive method, such as one focussed only on key immunogenic positions. Fourth, Distinctiveness can be sensitive to sequence alignment parameters. Complementary analyses that are independent of sequence alignments are warranted to overcome this shortcoming (44). Finally, Distinctiveness does not take into account amino acid similarities in the sequence alignments or the recency of the SARS-CoV-2 sequences used to build the alignment. Future work should account for amino acid similarities using substitution matrices (45) and incorporate the time of sequencing as parameters in computing the Distinctiveness scores.
In conclusion, we highlight that Distinctiveness more holistically captures the ongoing combat between viral evolution and host immunity, wherein lineages which are most distinctive from both the ancestral strain (the basis for all authorized COVID-19 vaccines) and VOCs (i.e. prior dominant strains against which natural immunity has developed) are the least likely to be neutralized by host immune responses. Distinctiveness can be considered as one important feature contributing to the epidemiological fitness of emerging SARS-CoV-2 variants, and thus, a salient factor to monitor as part of the global pandemic preparedness efforts.

Quantification of number of distinct positional amino acids for prevalent SARS-CoV-2 lineages
Individual substitutions, insertions, and deletions for each aligned SARS-CoV-2 protein sequence along with the corresponding PANGO designation were obtained from the GISAID (https://ww w.gisaid.org) database, on 2022 May 3. We considered only sequences labeled as "complete" and "high coverage" from the GI-SAID data, and collected from 28 top sequencing countries (Supplementary Table S1), this resulted in a total of 4,926,906 sequences. For the original, Wuhan strain and the five VOCs (Alpha, Beta, Gamma, Delta, and Omicron), the PANGO classification was obtained from the CDC website (https://www.cdc.gov/coronaviru s/2019-ncov/variants/variant-classifications.html).

Calculation of sequence Distinctiveness
For a given sequence, Distinctiveness within a geographical region of interest (i.e. a country) is defined as the average distances at the amino-acid level between that sequence and all sequences that were collected at least one calendar day before that sequence (limited by the time-resolution of the data). Specifically, for a sequence, s, it's Distinctiveness, D(s), is calculated using the following formula where N p is the number of prior sequences, s' is one specific prior sequence, the inner sum is over all pairwise aligned amino acid positions, and δ(s(p)-s`(p)) evaluates to 1 if sequence s and s`have the same amino-acid identity (one of 20 amino acids, a deletion, or a specific insertion) at position p and 0 otherwise. Positions of amino acids are determined relative to the Wuhan-Hu-1 reference, and insertions were treated as a single modification at the site of insertion. In cases where a nonsense mutation occurred, resulting in an early stop codon, mutations that followed this stop codon were not considered.

Calculation of sequence mutational load
The mutational load was calculated as the number of mutations away from the ancestral Wuhan-Hu-1 sequence. Similar to in the Distinctiveness calculation, insertions were counted as a single mutation. In cases where a nonsense mutation occurred, resulting in an early stop codon, mutations that followed this stop codon were not considered.

Calculating local prevalence of VOCs
The local prevalence of a SARS-CoV-2 variant, as reported in Figure 2 was calculated as the percentage of SARS-CoV-2 sequences in GISAID that were assigned to a lineage comprising that variant, during specific time windows and in specific countries.

Correlating the Distinctiveness and changes in future prevalence of SARS-CoV-2 lineages
We correlated the average Distinctiveness of sequences in a set during a 28-d window to the change in prevalence of the corresponding set, defined as prevalence (t + 56 to t + 84)-prevalence (t to t + 28), where t denotes time. For the analysis in Figure 3a Figure S4 and yields similar conclusions as those described in the main text. ROC-curves were generated from these data using Scikit-learn, using binary labels based on a minimum 20 percentage point increase in lineage prevalence for a country/time datapoint. Resulting AUC and threshold values, maximizing the sum of Sensitivity and Specificity, were found to be robust with respect to the cut-off used for labeling the data based on the percentage point increase (Supplementary Figure S5). We used bootstrap resampling (10,000 samples) of the underlying data points (scatter points in Figure 3a) to estimate 95% CIs on the resulting AUC values.

Labeling of neutralizing antibody epitope sites on the Spike protein
We have abstracted antibody epitope data for Therapeutic antibodies, as tracked by NCATS (https://opendata.ncats.nih.gov /covid19/), as well as Neutralizing antibodies, typically isolated from convalescent patient sera, as encountered in the Protein Data Bank (46). We define an antibody epitope as all residues in the antigen protein that have heavy (non-hydrogen) atoms at a distance of 4Å or less to heavy atoms of the bound antibody. When a structure of an antigen-antibody complex contains multiple instances of the interaction, such as in the case of a Spike protein trimer, and/or when several structures of the same antigen-antibody complex are available, we aggregate the binding data into a single epitope definition. We have also collected data for neutralizing antibodies as listed in Supplementary Data files provided by the Bloom and Xie labs (47,48), who have reported the results of single-point mutations that affect binding affinities (https://media.githubusercontent.com/media/jbloomlab/SARS 2_RBD_Ab_escape_maps/main/processed_data/escape_data.csv). We have listed residues whose mutations were found to have a nontrivial effect on binding activity for a given antibody (site total escape of 0.1 or higher). These are not necessarily close in 3D structure. As structures of those antibodies with bound antigen become available, we do find good agreement, in general, and we amend the epitope definition with that derived from the 3D structure data. In a few instances, structure-derived epitopes were slightly extended based on the characterization of the epitope by the structure's authors, and may include interactions slightly beyond the 4Å cutoff that we have employed.
Specifically, the following positions in the Spike protein were labeled neutralizing antibody epitope sites: 13

Supplementary material
Supplementary material is available at PNAS Nexus online.

Funding
This study was self-funded by nference. No external funding was received for this study. The work of A.M.B., J.W., R.C., and J.B. was supported by the National Center for Biotechnology Information of the National Library of Medicine, National Institutes of Health.

Authors' contributions
M.N. and K.M. designed research, performed research, analyzed data, and wrote the paper; P.L. designed research, analyzed data, and wrote the paper; A.M.-B. and J.W. wrote the paper; R.C. contributed new reagents/analytic tools, and wrote the paper; J.R.B. wrote the paper; A.J.V. designed research; analyzed data, and wrote the paper; and V.S. designed research, contributed new reagents/analytic tools, and wrote the paper.