key: cord-0429534-whtb5psq
authors: Bernasconi, Anna; Mari, Lorenzo; Casagrandi, Renato; Ceri, Stefano
title: Analysis of amino acid change dynamics reveals SARS-CoV-2 variant emergence
date: 2021-07-13
journal: bioRxiv
DOI: 10.1101/2021.07.12.452076
sha: be1930e73794cdac9eb473fb23b82116c6bf9b5a
doc_id: 429534
cord_uid: whtb5psq

Since its emergence in late 2019, the diffusion of SARS-CoV-2 is associated with the evolution of its viral genome1,2. The co-occurrence of specific amino acid changes, collectively named ‘virus variant’, requires scrutiny (as variants may hugely impact the agent’s transmission, pathogenesis, or antigenicity); variant evolution is studied using phylogenetics3–6. Yet, never has this problem been tackled by digging into data with ad hoc analysis techniques. Here we show that the emergence of variants can in fact be traced through data-driven methods, further capitalizing on the value of large collections of SARS-CoV-2 sequences. For all countries with sufficient data, we compute weekly counts of amino acid changes, unveil time-varying clusters of changes with similar – rapidly growing – dynamics, and then follow their evolution. Our method succeeds in timely associating clusters to variants of interest/concern, provided their change composition is well characterized. This allows us to detect variants’ emergence, rise, peak, and eventual decline under competitive pressure of another variant. Our early warning system, exclusively relying on deposited sequences, shows the power of big data in this context, and concurs to calling for the wide spreading of public SARS-CoV-2 genome sequencing for improved surveillance and control of the COVID-19 pandemic.

over 64 weeks, whereas South Africa has the smallest (84 changes, 56 weeks). We observed that, on emergence, the prevalence of changes that will eventually characterize a variant typically show an exponential-like temporal growth 30 ; one meaningful example is shown in Extended Data Fig. 1 .

We harnessed the peculiarity of these temporal trends to scout emerging variants. To that end, we applied standard time-series clustering techniques to group together changes showing similar prevalence behavior over a one-month-long period. We regarded as warnings of possible variant emergence those clusters characterized by (i) a positive trend in the prevalence time-series of the constituting amino acid changes, and (ii) being sufficiently different from clusters that caused previous warnings (to avoid extracting twice a set describing the same candidate variant). To assess whether these clusters successfully match known variants (named 'hit'), we tested their similarity against the characterizing changes of the most widespread SARS-CoV-2 lineages, as defined by Pangolin, collected in the 'lineage dictionary' (see Extended Data Tables 1 and 2). We also used community analysis to assess whether and how the emergence of new variants can lead to a reorganization in the dynamics of amino acid changes, as observed through their co-occurrence within clusters.

Our approach must not be considered an alternative to phylogenetic analysis, which takes into account the entire evolutionary history of the viral genome. Other methods have previously complemented phylogenesis, namely by describing typical SARS-CoV-2 mutational profiles across different countries and regions 31,32 , proposing statistical indicators for location-based mutation evolution 33 , and observing changes that become recurrently prevalent in different locations, thus suggesting selective advantages 1 . Time-based analyses have been considered for phenetic clustering of prevalent SARS-CoV-2 mutations over time 34, 35 , trend detection in SARS-CoV-2 short nucleotide sequences 36 , and single amino acid changes 37 .

To show our methods in action, Fig. 1 displays the case of Japan, featuring a matrix of 594 amino acid changes observed over 63 weeks. We identify 132 clusters, ten of which are included in the early warning system, tracking the emergence of possible variants. Four of these clusters revealed to be early hits, i.e. showing strong cluster-lineage similarity at the time of clustering. These hits occur on Jun-wk4-20, and March-wk3-21 (two instances), respectively associated with variants JP2, JP1, and JP3/Alpha. We use JP1, JP2, and JP3 for lineages B. 1.1.214, B.1.1.284, and R.1 since, according to Pango lineages reports 38 , they originated in Japan. Each hit is qualified by the growth of the change prevalences in the identified cluster (left panels a-c) and by the similarity between cluster composition and the lineage dictionaries (central panels a-c). Partial overlaps of identified clusters with more than one dictionary are due to overlaps between dictionaries. We find that the emergence of a new variant is typically associated with a drastic reorganization of the within-cluster co-occurrence patterns of amino acid changes (right panels a-c), which may suggest profound effects of variant emergence. Panels d-g of Fig. 1 show the temporal dynamics of the changes associated with each identified variant. Strong cluster-dictionary matches, illustrated by circle size and color, are found for all variants, spanning several months after the moments when warnings are issued. By contrasting variant dynamics with recorded cases (shaded bars), we observe that JP2, JP1, and Alpha prevalences peaked ahead of the first, second, and third COVID-19 waves in Japan. By contrast, JP3 shows a prevalence peak when reported cases were at a minimum. Taken all together, Fig. 1 constitutes a syntactic fingerprint of how variants have evolved in Japan. In addition to these four hits, our method identified other six candidates that could be rapidly dismissed, because clusters identified in subsequent weeks showed low similarity with the original candidate. Four of them were discarded two weeks after the warning, the other two within five weeks. Extended Data Fig. 2 shows all ten warnings.

Variant emergence around the world Fig. 2 illustrates how our method reconstructs the emergence of the main WHO-named variants at their country of origin. The temporal pattern of the Alpha variant is remarkably neat, exhibiting an initial phase of exponential growth, and well correlating with the number of reported COVID-19 cases. After a tenweeks plateau with prevalence values close to 100%, a five-weeks sharp decline occurred, associated with the concurrent growth of the Delta variant (see Extended Data Fig. 3 ). Note that some variants never reached a high prevalence, e.g. Epsilon and Iota. Variant dynamics also differed in terms of peak timing. For instance, variant Zeta well matched the temporal evolution of the second COVID-19 wave in the country, and Gamma shortly anticipated the third large wave of infections and became locally dominant. Similarly, variants Kappa, Delta, and Beta, grew quite significantly before the second wave of cases in their respective countries of origin. Note that some weaker matches are found also before the relevant warnings, essentially because some amino acid changes included in the lineage dictionaries did indeed emerge individually well before the diffusion of variants. Also, note that the WHO-named variants Eta, Theta, and Lambda do not appear in this figure as their countries of origin did not meet our criteria for minimum data availability and were not included in the analysis. Fig. 3 shows the temporal patterns of notable variants as identified by our method in Europe and in the US. For all European countries with at least 16k deposited sequences, panel a displays the growth curves of changes associated with the Alpha variant at the time when the relevant warnings were issued by our method. The spread of Alpha outside the UK was delayed by five to seven weeks; warnings concentrate on three consecutive weeks in late January (country shading). Interestingly, the growth curves of the change prevalences associated with the Alpha variant outside the UK are more noisy than inside the UK, likely due to availability of fewer sequences and spatial heterogeneities in sequence deposition.

Panel b of Fig. 3 shows instead the dynamics of the most widespread US-born SARS-CoV-2 variants in all US states where at least 5k sequences were deposited and their interplay with WHO-named variants.

Our clustering and warning methods partitioned these states into three distinct, geographically quite compact, groups. A central/eastern group of states (blue) displays a baseline temporal pattern where common US variants first rose until Jan-wk1-2021, then started declining under pressure of variant Alpha. In the Western group of states (orange), the Californian-born variant Epsilon appeared in-between the US2 variant and Alpha, whereas in the north-eastern group (green) the New-York-born Iota appeared after the rise of Alpha, none of them reaching the same prevalence as Alpha. Although quite general, the patterns above described in groups may show country-dependent peculiarities. In Louisiana, for example, the hit for Epsilon followed the one for Alpha, and we did not issue warnings for Alpha in New York, nor for US1/US2 in New Jersey and Connecticut.

We have presented a novel methodology to scout SARS-CoV-2 variants entirely based on the analysis of data from the GISAID sequence repository; specifically, we have shown that emerging lineages can be detected via cluster analyses of the prevalence time-series of amino acid changes. Our scouting method is not in competition with phylogenetic analysis, which is fundamental for properly defining variants as lineages through monitoring of the evolution of viral sequences. Nevertheless, the proposed big-data-driven approach proved to be quite effective. First, the major WHO-named variants of SARS-CoV-2 were well and timely captured in their country of origin, an indication that the early dynamics of change prevalence is a distinguishing property of emerging variants. All hits were fast (requiring no more than five weeks, the length of the time window used for cluster analysis), except the hits of Beta (five plus two weeks) and Iota (five plus one weeks)-see Table 1 . Second, our plots reveal that the method can effectively trace and visualize the variants' natural evolution. Variants emerge, reach their maximum prevalence, and eventually start a declining phase, typically when a competing variant emerges. Our method also issues warnings that are not found to match major lineages down the road. We do not see this outcome as an intrinsic limitation of our approach, though. In fact, even when not selected as hits, warnings are informative of the dynamics underlying variant emergence and inter-variant competition; in addition, similarity analysis can clarify, typically in a matter of a few weeks at most, whether an emerging variant is gaining traction.

Relying exclusively on big data for variant scouting is hampered by three problems: (i) consistency of sampling, (ii) delay of deposition, and (iii) biases in sampling. As for (i), the ratio between the number of sequences and reported COVID-19 cases widely varies among countries and drops from 77% in Iceland (clearly facilitated by small number of cases) to below 0.1% for many countries, also including some large ones like India and Brazil. Among the countries with lots of cases, the UK stands with an exceptionally high ratio exceeding 9%. Indeed, US and UK have contributed the largest number of sequences worldwide (517k and 425k, respectively). Extended Data Table 3 shows these statistics for all countries contributing to GISAID with more than 1,000 sequences, whereas Extended Data Table 4 shows statistics for the 50 US states. Regarding (ii), as an example, in the UK the average delay between collection and deposition amounts to 24 days. This delay tended to reduce as the pandemic unfolded, from 38 days in 2020 to just 16 days in 2021. Iceland is again striking the best performance, with 11 days of average delay in 2021. Concerning (iii), heterogeneities in surveillance may introduce biases that should be considered, for instance when sampling is concentrated within regions where variants have newly appeared, thereby causing possible over-estimation of the observed prevalence of variants within a country. Another aspect interfering with a correct estimation of variants' prevalence is the concurrent rollout of COVID-19 vaccines, particularly in case of variants possibly endowed with partial immune escape potential, as in the recent case of Delta variant 39 .

Being aware of these obstacles, the potential of developing and using big-data-driven approaches to keep track of variant emergence may be proved by comparing the timings of warnings issued by our method for the eight major WHO-named variants (as of June 3rd, 2021) with the dates of their initial communication in the media or in relevant scientific outlets (see again Table 1 and Methods). In five cases (Alpha, Beta, Epsilon, Zeta, Kappa), our data-driven warnings anticipated the press notifications, in some cases (Beta, Epsilon, Kappa) with a lead time of several weeks. In other two cases (Delta, Iota), the warning would have been issued around the same time as the first communications. Only in the case of the Gamma variant our warning was late (by a couple of weeks) compared to its first communication. It must be remarked that Brazil, where both Gamma and Zeta variants originated, was sequencing a low number of samples at that time.

In conclusion, our warnings were issued at times which favourably compare with most press/research announcements, in some cases even anticipating them. The above figures illustrating variant hits and dynamics collectively provide an effective fingerprint of the location-specific temporal evolution of variants, and well support comparative displays in time and space. Our exercise thus strongly corroborates the urgent call, raising from the entire scientific community, for faster and wider public publishing of viral sequences 40 . ; values J cd < 0.1 are not shown. Right panels: Community detection applied to the within-cluster change co-occurrence matrix evaluated over the period of time between the emergence of two subsequent variants. The resulting interaction network is plotted using a force-directed layout 41 , with edge lengths inversely proportional to link weights. Colors indicate nodes belonging to communities that strongly match known lineages, other communities are in gray-scale shades. d-g, Temporal dynamics of the four identified variants. In each scatter plot, the left y-axis value represents the average prevalence of the changes in the cluster-dictionary intersection, circle size represents the cluster-dictionary similarity (J cd < 0.1 not shown), while circle color is proportional to the log-ratio between the number of changes in the cluster and in the dictionary. Warnings are marked with a labeled arrow. In the background, the number of reported COVID-19 infections (thousand cases, right axis) in Japan 42 . location, and the list of amino acid changes. All records not reporting the day (but only year, or year and month) in the collection date field were discarded, resulting in a dataset of 1,763,923 records (about 97% of the initial size). To denote amino acid changes, we adopt the common notation used by GISAID, representing them as the concatenation of the following elements: the protein acronym, the reference amino acid residue, the position (using the coordinate system of the single protein), and the alternative amino acid residue exhibited by the specific mutated sequence (in case of substitutions) or a dash (in case of deletions). Proteins include, for example, Spike, Envelope (E), Membrane (M), Nucleocapside (N), as well as non-structural proteins forming the ORF1ab polyprotein (NSP1, NSP2, . . . , NSP16). Changes can occur at any point within a protein; they are evaluated with respect to the reference sequence WIV04 (used by GISAID 47 ).

Original collection dates were temporally binned into four approximately weekly periods per month, namely days 1-7 (shortcut as wk1), 8-15 (wk2) , 16-23 (wk3), and 24-end of month (wk4). We verified that the small differences in bin sizes among the four periods do not have any impact on the results obtained using our data mining methods, as they rely on change prevalence (relative count of sequences) rather than change abundance (absolute count). We transformed the GISAID data into a table containing tuples of the following kind: amino acid change, country of collection, week of collection, absolute count of sequences holding that change, and total collected sequences in the same country-week. As locations of interest, we selected all European countries and US states for which at least 16k and 5k sequences were available, respectively, and a small number of other countries where variants of concern/interest have originated according to the World Health Organization 17 . For each of those countries, we prepared a matrix where: each row r represents an amino acid change that has been observed in at least five sequences and for at least three weeks, each column c represents an observation week, each cell r, c contains the number of sequences exhibiting the amino acid change r observed in week c. The full account of matrix sizes for the selected countries is shown in Extended Data Table 5 . We also extracted from the GISAID dataset the total number of sequences collected at given locations for each of the considered weekly periods. Data extraction and aggregation was performed using 

We extracted all lineages that are expressed in at least 5k sequences worldwide and a small number of lineages that deserved attention by the news and/or other online sources, provided they appeared in at least 1k sequences worldwide. The assignment of amino acid changes to lineages is not uniquely/uniformly defined by all sources; following the choice made by Mullen et al. 45 , we computed the set of characterizing changes, pragmatically defined as those that appear in at least 75% of the lineage sequences in GISAID, for all retained lineages. This rule produces partially overlapping lists of amino acid changes, with lengths ranging from five to 24, which we regard as our 'lineage dictionary' (Extended Data Table 2 ). To improve readability, the naming convention starts with the Pangolin lineage 14 , followed by the WHO name 17 (using Greek alphabet) when available. In selected cases, we add a geographical characterization, based on the country of origin (https://cov-lineages.org/ 38 ), resulting into US1/US2 and JP1/JP2/JP3 aliases. Only exceptionally we merged multiple lineages into one label, namely when lineages share the same set of characterizing amino acid changes. Extended Data Table 1 contains a comprehensive list of lineages considered in our study.

We found that few changes (globally six, four worldwide -S_D614G, NSP12_P323L, N_R203K, and N_G204R -and two specific to North-America-NS3_Q57H and NSP2_T85I, out of a total of 78,650 recorded changes) are disproportionately frequent in the dataset of each continent, being present in more than 60% of GISAID sequences and populating the dictionary of more than 30% of the most frequent lineages (i.e. those represented by at least 80 sequences). To avoid confounding effects in data analysis, we removed these changes from both the aggregated data matrices and the dictionary entries.

To assess whether the temporal patterns of changes' prevalence occur in a coherent manner in the dataset, as expected in the case of changes belonging to closely related variants, we apply a relatively simple time-series clustering algorithm. All the following analyses have been performed using MATLAB R2021a. At each time t, we retain the n(t) time-series of change prevalence (current ratio between change counts and total counts of changes) that are observed continuously over a time interval of four weeks prior to t ([t − w + 1 · · · t], with w = 5) and partition them via k-medoids clustering 48 (PAM algorithm, kmedoids function in MATLAB), with pairwise distances between time-series being evaluated via dynamic time warping 49 (dtw in MATLAB). The optimal value of k is exhaustively searched over the range [1 . . . min(20, n(t 

selected as the one that maximizes the average silhouette score 50 evaluated over the entire set of n(t) change prevalence time-series.

To identify clusters of changes characterized by an increasing trend in their prevalence time-series, as expected in the case of emerging variants, we use the non-parametric Kendall's τ B statistic 51 , as implemented in the ktaub MATLAB package 52 , namely to evaluate the ordinal association between change prevalence and sampling time. More precisely, a cluster of changes is considered to deserve further attention as a candidate for variant emergence if its average prevalence time-series shows a positive trend (τ B > 0 at significance level α = 0.05). Among these candidate clusters, we are particularly interested in those that are sufficiently different from previously observed trending clusters, because they could signal the emergence of a new virus variant. Indeed, the partitioning of change time-series is performed independently at each time step, thus the clusters identified at a given time are in principle not related to those identified at any previous step. However, because of the temporal autocorrelation of change prevalence dynamics, similarities are expected (and found) to exist in the composition of clusters identified at subsequent time steps. To quantify similarity for each pair of clusters, we use the classical Jaccard index 53 , defined as the ratio J cc between the cardinality of the intersection and the cardinality of the union of the sets of changes constituting the two clusters. Among the candidate clusters with a positive trend in their prevalence time-series, only those that are sufficiently different from any previous candidates (J cc < 0.5) are selected to form an early warning system for monitoring the possible emergence of new variants.

We use the similarity between the composition of early warning clusters and the lineage dictionaries to assign a posteriori the observed change time-series to known lineages, thereby building a ground truth for assessing the performance of data analysis. Specifically, we use again Jaccard similarity index, this time applied to cluster vs. dictionary change composition (J cd ). A threshold J cd > 0.5 is used to identify clusters for which there is a close compositional match to a known lineage. If a match exists, we conclude that the method has identified the lineage (thus, the underlying variant) at the same time as when clustering occurred (early hit). If not, we keep monitoring the cluster-dictionary similarity J cd in the following time steps, each time using the cluster with the highest compositional similarity J cc with that of the previous step, provided that the inter-step similarity remains sufficiently strong (J cc > 0.5), until the threshold J cd = 0.5 is possibly exceeded (late hit). If the inter-step similarity test fails without producing a late hit, the candidate cluster causing the early warning is discarded. Rapid convergence towards a decision between late hit or candidate discarding constitutes a desirable property of our method.

The study of the changes associated with the successfully identified variants can be deepened by analyzing their temporal dynamics along with the country-specific epidemiological patterns of the pandemics.

Specifically, we match the dictionary composition of each identified variant to cluster composition over time and evaluate: i) the average prevalence at the time of clustering of the changes in the intersection between cluster composition and lineage dictionary; ii) the cluster-dictionary similarity J cd ; and iii) the log-ratio between cluster and dictionary cardinalities. This methodology is applied to the whole dataset, i.e. also prior to the emergence of variants or after their possible disappearance; in fact, some changes belonging to a lineage dictionary may be individually present in the population before/after the diffusion of the relevant variant.

We collected information about the first Institutional Communications (ICs), Research Communications (RCs) captured at their first version, and Published Papers (PPs) about the most important WHO-named variants to assess whether the points in time when warnings are issued by our method compare well with the points in time when variants actually became known. The following references provide the dates indicated in Table 1 : 

The points in time when lineage dictionaries are for the first time successfully (J cd > 0.5) matched to cluster composition are used to partition temporally the dataset of change prevalence. For each resulting time window, we explore the within-cluster co-occurrence dynamics of the different changes by building a weighted, undirected graph where the nodes are the different changes and the edges represent pairwise interactions between changes; edge labels are evaluated as the number of clusters in which two changes occur robust co-occurrence patterns within the clusters identified over the time window of interest, we apply the Louvain method of community detection 66 using the MATLAB function GCModulMax1 (Community Detection Toolbox 67 ). The change compositions of the resulting communities are compared with the lineage dictionaries by using again Jaccard index.

All data matrices used in this study, one for each country or US state, as well as the lineage dictionary, have been deposited in Zenodo at https://doi.org/10.5281/zenodo.5090169. Original sequence metadata and amino acid changes are publicly accessible through the GISAID platform.

The Zenodo repository https://doi.org/10.5281/zenodo.5090169 includes custom scripts that can be used to reproduce the results presented in this study. Extended Data Figure 2 | Data-driven warning of possible candidate variant emergence issued in Japan. a, The nine time points when warnings were raised. Note that Mar-wk1-21 exhibits two distinct warnings for clusters JP_093 and JP_98. b, Cluster-lineage Jaccard similarity for the ten warnings. The four hits (similarity above the threshold J cd = 0.5) are reported in Fig. 1 in the main text. Strong cluster-lineage similarity is already observed at the time when clustering is performed, which qualifies these as early hits. In this case, candidates were discarded within five weeks, as no late hits (i.e. strong similarity observed with some delay after clustering/warning between a cluster with an amino acid change composition that can be traced back to the original warning and a lineage dictionary) were recorded. Figure 3 | Temporal dynamics of the two variants identified in UK. We show the decline of Alpha variant occurring together with the growth of the Delta variant. Details as in panels d-g of Fig. 1 Extended Data Figure 4 | Data extraction and aggregation steps. GISAID original data comprise records of sequence accession ID, collection date, submission date, Pangolin lineage, collection location, and list of amino acid changes. We compute a table of tuples with: amino acid change, country of collection, week of collection, absolute count of sequences holding that change, and total collected sequences in the same country-week. For selected countries (the example shows Japan), we build a matrix where rows are amino acid changes observed in at least five sequences for at least three weeks, columns are observation weeks, and cells contain the number of sequences collected in the given week exhibiting the specific amino acid change. A separate vector holds the total number of sequences collected in the location under study for each of the considered weekly periods. Note that the last weeks of May 2021 hold very few data, thereby producing diluted results of our clustering and warning method. Countries contributing to the GISAID database with more than 1,000 sequences; COVID-19 cases 42 ; percentage of sequences over cases; average delay (in days) between sample collection and sequence submission dates.

Tracking changes in SARS-CoV-2 Spike: evidence that D614G increases infectivity of the COVID-19 virus

One year of SARS-CoV-2 evolution

Phylogenomics reveals viral sources, transmission, and potential superinfection in early-stage COVID-19 patients in Ontario

COVID-19 in Amazonas, Brazil, was driven by the persistence of endemic lineages and P.1 emergence

Genomic epidemiology of SARS-CoV-2 reveals multiple lineages and early spread of SARS-CoV-2 infections in Lombardy, Italy

Spread of a SARS-CoV-2 variant through Europe in the summer of 2020

Global initiative on sharing all influenza data-from vision to reality

Genetic variants of SARS-CoV-2-what do they mean?

SARS-CoV-2 samples may escape detection because of a single point mutation in the N gene

Mutations on COVID-19 diagnostic targets

Efficacy of the ChAdOx1 nCoV-19 Covid-19 vaccine against the B.1.351 variant

Sensitivity of infectious SARS-CoV-2 B.1.1.7 and B.1.351 variants to neutralizing antibodies

Multiple SARS-CoV-2 variants escape neutralization by vaccine-induced humoral immunity

Emergence of a novel SARS-CoV-2 strain in Southern California, USA

Genomic characterization of a novel SARS-CoV-2 lineage from Rio de Janeiro

A novel SARS-CoV-2 variant of concern, B.1.526, identified in New York

Ministry of Health and Family Welfare, India. Genome Sequencing by INSACOG shows variants of concern and a Novel variant in India

Last accessed

The time varying network of urban space uses in Milan

Fast unfolding of communities in large networks

Retrieved

Extended Data Figure 1 | Counts of protein Spike amino acid changes in the UK in different bi-weekly periods

A blue background is used in the positions notably belonging to the Alpha variant

higher black bars in these positions from top to bottom indicate an increase in the frequency of these amino acid changes. B.1.1.214(JP1) N_M234I

B.1.1.284(JP2) N_P151L, NSP12_A423V

B.1.1.519 NSP3_P141S, NSP4_T492I, NSP6_I49V, NSP9_T35I, Spike_P681H

B.1.1.7(Alpha)

B.1.177.86 N_A220V, NS7b_E39*

B.1.2(US1) N_P199L, N_P67S, NS3_G172V, NS3_Q57H, NS8_S24L, NSP14_N129D, NSP16_R216C

B.1.221 N_P199L, NS3_G172R, NS3_Q38R, NS3_V202L

B.1.258.17 NS3_Q185H, NS8_E64*, NSP12_V720I, NSP13_A598S, NSP13_H290Y, NSP13_P53L, NSP14_E453D, NSP14_P46L, NSP3_I1683T

B.1.427(Epsilon) NS3_Q57H, NSP13_D260Y, NSP13_P53L, NSP2_T85I, NSP4_S395T, N_T205I, Spike_L452R

B.1.429(Epsilon) NS3_Q57H, NSP13_D260Y, NSP2_T85I, NSP9_I65V, N_T205I, Spike_L452R

Spike_Y144-B.1.526(Iota) N_M234I, N_P199L, NS3_P42L, NS3_Q57H, NS8_T11I, NSP13_Q88H, NSP2_T85I, NSP4_L438P

Spike_Y144-B.1.526.2(Iota) N_P13L, N_S202R, NS3_P42L, NS3_Q57H

B.1.596(US2) N_P199L, N_P67S, NS3_G172V, NS3_Q57H, NS8_S24L, NSP14_N129D

B.1.617.1(Kappa) N_D377Y, N_R203M, NS3_S26L, NS7a_V82A, NSP13_M429I, NSP15_K259R, NSP3_T749I, NSP6_T77A, Spike_E484Q, Spike_L452R

B.1.617.2(Delta)

P.1|P.1.2(Gamma)

Zeta) N_A119S, N_M234I, NSP5_L205V, NSP7_L71F

R.1(JP3) M_F28L, N_Q418H, N_S187L, NSP13_G439R, NSP14_P412H, Spike_E484K

The dictionary contains amino acid changes occurring in at least 75% of a lineage's sequences. Worldwide confounding (too frequent) changes are already removed, North-America-specific ones are not

Extended Data Table 4 | Statistics on SARS-CoV-2 sequencing in all US states

We gratefully acknowledge all data contributors, i.e. the Authors and their Originating Laboratories responsible for obtaining the specimens, and their Submitting Laboratories that generated the genetic sequence and metadata and shared via the GISAID Initiative the data on which this research is based. A.B. and S.C. acknowledge the support provided by the ERC Advanced Grant 693174 "Data-Driven Genomic Computing (GeCo)".

A.B. and S.C. proposed the study, all authors jointly conceptualized and designed it; A.B. and S.C. designed the data extraction methods (implemented by A.B.); R.C. and L.M. designed the data clustering and analysis methods (implemented by L.M.); A.B. mined and curated published SARS-CoV-2 variant data; all authors wrote and revised the manuscript.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to A.B. (anna.bernasconi@polimi.it